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ABSTRACT 

This  report  describes  the  target  detection  performance  of  an  airborne  surface  surveillance  radar  in 
the  presence  of  medium  grazing  angle  sea-clutter.  In  the  absence  of  frequency  agility,  the  temporal 
correlation  of  the  sea-clutter  can  be  significant  and  if  it  is  not  accounted  for  in  the  radar  model, 
the  required  signal  to  interference  ratio  for  a  given  probability  of  detection,  P&,  will  be  incorrect 
by  several  dB,  resulting  in  over-estimated  performance.  This  report  describes  a  robust  method  for 
calculating  the  Pd  for  both  K  and  Pareto  compound  sea-clutter  distributions.  Empirical  models  of 
the  amplitude  distribution  and  the  speckle  correlation  are  used  to  determine  the  expected  detection 
performance  given  different  collection  geometries  and  environmental  conditions  with  the  output 
used  to  determine  the  minimum  detectable  target  radar  cross  section  in  a  detection  scenario. 
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Radar  Detection  Performance  in  Medium  Grazing  Angle 

X-band  Sea-clutter 


Executive  Summary 


This  report  summarises  the  target  detection  performance  of  an  X-band  surface  surveillance  radar 
operating  at  medium  grazing  angles  in  the  maritime  environment.  A  number  of  different 
detection  scenarios  are  explored  with  variations  in  the  collection  geometry,  sea-state  and 
polarisation.  Two  sea-clutter  distributions  are  considered:  the  popular  K-distribution  and  the 
Pareto  distribution  which  has  recently  been  shown  to  be  a  better  match  to  the  distribution  tail. 
Consequently,  detection  performance  results  using  the  Pareto  distribution  have  been  chosen  as 
‘truth’  in  the  report. 

To  test  these  scenarios,  a  number  of  parameter  models  are  required  to  relate  the  radar  and 
environmental  characteristics  to  the  mean  backscatter  and  distribution  shape  values.  Appropriate 
models  are  therefore  presented  along  with  a  full  mathematical  description  of  the  different  target 
fluctuation  models.  A  major  contribution  of  the  report  is  a  full  derivation  of  the  Ward,  Tough  and 
Watts  detection  algorithm  for  correlated  clutter  with  a  modification  to  make  the  technique  robust 
for  high  clutter  to  noise  ratios. 

The  sea-clutter  and  target  detection  models  are  then  used  to  analyse  the  probability  of  false  alarm, 
Pfa,  probability  of  detection,  I\\  and  the  minimum  detectable  radar  cross  section  (RCS)  for  a  given 
set  of  radar  parameters.  The  Pfa  results  show  a  big  difference  in  the  horizontal  polarisations  (HH 
and  HV)  with  the  K  distribution  underestimating  the  threshold  by  up  to  5  dB,  while  there  is  only 
a  minor  difference  in  the  vertical  channel  (VV).  The  effect  of  non-coherently  integrating  pulses  is 
to  reduce  the  threshold. 

A  comparison  of  the  target  fluctuation  models  reveals  that  the  Pd  curve  is  steeper  for  the  chi- 
squared  distribution  model  than  the  exponential  model.  Increasing  the  target  correlation  reduces 
the  slope  for  both  models,  with  the  effect  being  more  pronounced  for  the  exponential  model.  Using 
a  distribution  model  also  allows  us  to  extrapolate  down  to  lower  Pfa  levels  than  could  be  used  in  a 
real  detection  scheme.  By  comparing  Pfa  levels,  it  is  found  that  for  the  HH  polarisation,  there  is 
a  larger  difference  in  the  signal  to  interference  ratio  between  the  K  and  Pareto  distributions.  This 
mismatch  also  got  larger  as  the  Pfa  reduced. 

Further  analysis  of  the  detection  performance  with  varying  azimuth,  grazing  and  sea-state  reveals 
that  the  choice  of  model  is  critical  to  accurately  determine  the  detection  performance.  There 
is  significant  mismatch  between  the  K  and  Pareto  distribution  results,  particularly  for  the  HH 
polarisation.  Other  large  discrepancies  are  observed  in  the  cross  wind  direction  and  at  lower 
grazing  angles.  The  mismatch  for  the  VV  polarisation  is  typically  very  small. 
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The  final  result  looked  at  the  minimum  detectable  target  RCS  given  a  set  of  radar  parameters.  For 
the  lower  grazing  angles  in  the  HH  and  F1V  channels,  there  was  a  large  difference  between  the  K 
and  Pareto  distributions,  which  reduced  as  the  grazing  angle  increased.  The  two  distributions  gave 
almost  identical  results  for  the  VV  channel.  In  terms  of  detectability,  the  Pareto  result  shows  that 
the  HV  polarisation  has  the  smallest  target  RCS,  then  VV  and  HH. 
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1  Introduction 

This  report  is  concerned  with  quantifying  the  detection  performance  of  an  X-band  surface 
surveillance  radar  in  the  maritime  environment.  The  results  are  based  on  data  from  DST  Group’s 
Ingara  radar  which  has  collected  12  days  of  fine  resolution  fully  polarimetric  sea-clutter  data  in  a 
circular  spotlight  mode  covering  all  azimuth  angles  and  focussing  on  the  medium  grazing  angle 
region,  15°  —  45°.  This  data  has  been  used  to  characterise  the  statistical  properties  of  the 
sea-clutter  in  order  to  accurately  quantify  target  detection  performance.  These  include  the  mean 
radar  backscatter  [Crisp  et  al.  2008],  amplitude  distribution  [Dong  2006,  Rosenberg,  Crisp  & 
Stacy  2010ft,  Weinberg  2011a,  Rosenberg  &  Bocquet  2015,  Rosenberg,  Watts  &  Bocquet  2014] 
Doppler  spectrum  [Rosenberg  2014]  and  sea  spikes  [Rosenberg  2013]. 

These  studies  then  led  to  a  number  of  empirical  models  which  take  into  account  different 
combinations  of  the  polarisation,  collection  geometry  and  environmental  conditions.  Relevant 
models  include  the  imaging  radar  systems  group  (IRSG)  mean  backscatter  model  [Crisp 
et  al.  2008],  the  short-time  temporal  correlation  [Rosenberg  2014]  and  the  K  distribution  shape 
parameter  [Rosenberg,  Crisp  &  Stacy  2009,  Crisp,  Rosenberg  &  Stacy  2015].  To  effectively 
model  the  effect  of  sea-spikes,  the  latter  model  has  also  been  modified  and  applied  to  the  Pareto 
distribution  shape  [Rosenberg  et  al.  2015].  These  empirical  models  are  described  in  Section  2  of 
the  report. 

Parametric  modelling  is  a  useful  tool  to  predict  a  radar’s  target  detection  performance  given  a 
statistical  description  of  the  environment.  Once  a  clutter  model  has  been  chosen  and  its 
parameters  determined,  a  threshold  can  be  determined  based  on  the  desired  probability  of  false 
alarm,  Pfa,  to  ensure  a  CFAR  (Constant  False  Alarm  Rate).  Any  backscatter  response  above  the 
threshold  is  then  declared  to  be  a  potential  target.  Using  this  threshold,  a  probability  of  detection 
Pd  can  be  calculated  for  a  given  target  fluctuation  model  to  determine  the  required  signal  to 
interference  (clutter  and  noise)  ratio  (SIR).  Subsets  of  the  empirical  parameter  models  have 
previously  been  used  to  analyse  target  detection  performance  in  K  and  KK  distributed  sea  clutter 
[Rosenberg  &  Crisp  2010a]  and  to  study  the  impact  of  temporal  correlation  with  different 
fluctuating  targets  [Rosenberg  2012,  Rosenberg  &  Bocquet  2013ft].  These  studies  found  that  if 
the  correct  distribution  was  not  used  to  model  the  tail  of  the  distribution  and  if  the  temporal 
correlation  was  not  accounted  for,  then  the  required  SIR  for  a  desired  Pd  can  be  over-estimated 
by  up  to  several  dB. 

The  main  contribution  of  this  report  to  provide  a  comprehensive  mathematical  description  of  the 
modelling  framework  and  to  extend  the  previous  results  using  the  newly  developed  Pareto  shape 
model.  Section  3  describes  the  modelling  framework  for  different  fluctuating  targets  using 
different  sea-clutter  distributions.  A  full  derivation  of  the  Ward,  Tough  &  Watts  [2013]  detection 
algorithm  for  correlated  clutter  and  /or  targets  is  then  presented  with  the  modification  [Rosenberg 
&  Bocquet  2013ft]  which  makes  the  technique  robust  for  high  clutter  to  noise  ratios  (CNR). 

The  final  Section  4  then  uses  the  empirical  models  to  quantify  the  detection  performance  for  a 
number  of  different  geometries  and  environmental  conditions.  Due  to  its  good  fit  to  the  Ingara 
sea  clutter  data  [Rosenberg  &  Bocquet  2013a],  the  Pareto  distribution  is  used  as  a  proxy  for  the 
data  with  the  K-distribution  results  showing  the  mismatch  that  would  be  expected  in  a  detection 
scenario.  These  results  are  then  fed  into  the  radar  range  equation  to  determine  the  minimum 
detectable  radar  cross  section  (RCS)  for  a  given  set  of  radar  parameters. 
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2  Sea-clutter  models 

To  comprehensively  characterise  target  detection  performance,  [Ward,  Tough  &  Watts  2013]  state 
that  knowledge  of  the  following  sea-clutter,  radar  processing  and  environmental  characteristics  is 
required: 

•  mean  backscatter 

•  amplitude  statistics 

•  temporal  correlation  or  spectrum 

•  spatial  characteristics  and  correlations 

•  discrete  clutter  spike  characteristics 

•  carrier  frequency 

•  polarisations 

•  spatial  resolution 

•  use  of  frequency  agility 

•  grazing  angle 

•  wind-speed  and  direction 

•  sea-swell  and  direction 

To  this  end,  a  number  of  experiments  were  conducted  using  the  Ingara  X-band  radar  to  collect  sea- 
clutter  covering  as  many  different  combinations  of  these  parameters  as  possible.  A  description  of 
the  trials  and  the  collected  data  is  presented  in  Section  2.1. 

Detection  performance  requires  an  accurate  representation  of  the  amplitude  or  intensity 
probability  distribution  function  (PDF).  The  most  commonly  used  PDF  model  for  sea-clutter  is 
the  K-distribution  as  it  captures  the  bulk  of  the  sea-clutter  PDF  very  well.  It  was  first  applied  to 
modelling  probability  distributions  of  sea-clutter  by  Jakeman  &  Pusey  [1976]  and  was  later  put 
into  a  Bayesian  or  compound  formulation  by  Ward  [1981].  This  allows  a  more  meaningful 
understanding  of  the  two  main  components  -  these  being  the  temporal  or  fast  varying  component 
and  the  slowly  varying  texture  component.  The  fast  component  is  commonly  known  as  speckle 
and  is  the  result  of  constructive  and  destructive  interference  effects  between  multiple  scatterers.  It 
is  typically  associated  with  the  small  local  wind-driven  ripples  (capillary  waves)  on  the  ocean 
surface.  The  texture  component  represents  changes  in  the  large  and  medium  scale  waves  which 
modulate  the  speckle.  In  terms  of  scattering  processes,  reflections  off  the  wind-driven  ripples  are 
known  as  Bragg  (or  resonant)  scattering,  while  the  non-Bragg  scattering  can  be  described  as  a 
combination  of  whitecaps  and  discrete  sea-spikes  [Ward,  Tough  &  Watts  2013].  It  is  these 
sea-spikes  which  cause  the  K-distribution  model  to  break  down,  particularly  at  finer  resolutions 
where  the  average  backscatter  reduces  and  the  sea-spikes  cause  higher  levels  of  backscatter 
relative  to  the  mean,  hence  extending  the  tail  of  the  distribution. 

In  order  to  better  model  this  effect,  the  KA  [Middleton  1999],  KK  [Dong  2006,  Rosenberg,  Crisp  & 
Stacy  20106],  Pareto  [Bailed,  Nehorai  &  Wang  2007,  Farshchian  &  Posner  2010,  Weinberg  201 16, 
Rosenberg  &  Bocquet  2015]  and  K+Rayleigh  [Lamont-Smith  2000,  Rosenberg,  Watts  &  Bocquet 
2014]  models  have  been  proposed.  The  KA  model  was  originally  proposed  by  Middleton  [1999] 
and  then  specifically  applied  to  spiky  sea  clutter  by  [Ward  &  Tough  2002,  Ward,  Tough  &  Watts 
2013],  A  comparison  of  this  model  and  the  KK  distribution  has  been  applied  to  the  Ingara  data 
[Dong  2006].  The  KK  distribution  offered  a  better  fit,  but  relied  on  fitting  up  to  five  parameters. 
A  number  of  optimal  and  sub-optimal  detectors  have  also  been  proposed,  but  they  are  difficult 
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to  implement  in  practice  [Weinberg  2012a].  This  led  to  the  Pareto  model  which  is  described 
by  two  parameters,  yet  can  reasonably  model  the  long  tails  present  in  the  sea-clutter  distribution. 
The  Pareto  distribution  also  has  significantly  simpler  optimal  and  sub-optimal  detectors  [Weinberg 
2012/?,  Weinberg  2014]. 

Clearly  the  PDF  must  account  for  the  different  elements  of  the  sea-clutter  as  well  as  the  thermal 
noise  in  the  radar.  Therefore,  the  first  element  required  for  modelling  is  the  CNR  which  is 
described  in  Section  2.2.  This  is  determined  using  the  thermal  noise  power  and  clutter  power 
(radar  range)  equations  [Ward,  Tough  &  Watts  2013]  with  the  IRSG  mean  backscatter  model 
[Crisp  et  al.  2008]  used  to  provide  variation  with  the  viewing  geometry,  sea  state  and  polarisation. 
The  K  and  Pareto  distributions  are  then  described  in  Section  2.3,  with  details  on  how  each  of  the 
shape  parameters  were  estimated.  The  newly  developed  shape  model  for  these  distributions 
[Rosenberg  et  al.  2015]  is  also  presented  in  this  section. 

To  improve  the  probability  of  detection,  a  radar  will  typically  average  returns  from  multiple  looks 
-  pulses  and/or  scans.  This  reduces  the  variance  of  the  random  interference  and  increases  the  mean 
separation  with  the  target  signal  component.  In  the  absence  of  frequency  agility,  there  may  be  little 
difference  in  the  radar  return  when  looking  at  the  same  clutter  patch  closely  in  time.  This  effect 
is  known  as  temporal  correlation  and  is  often  characterised  from  its  Fourier  transform  or  in  this 
context,  its  radar  Doppler  spectrum.  There  is  a  large  body  of  literature  focussed  on  understanding 
the  radar  Doppler  spectrum  [Long  2001,  Ward,  Tough  &  Watts  2013,  Plant  1997].  From  this 
analysis  came  two  significant  models  describing  the  line  shape  of  the  mean  Doppler  spectrum 
from  low  grazing  angles  [Lee  et  al.  1995,  Walker  2001].  The  former  uses  a  Gaussian  component 
to  model  the  Bragg  scatterers  and  a  Lorentzian  and/or  Voigtian  for  the  non-Bragg  components.  A 
simpler  model  was  later  presented  by  Walker  [2001],  who  uses  a  combination  of  three  Gaussians 
to  describe  the  Bragg,  whitecap  and  sea-spike  components  of  the  Doppler  spectrum.  At  medium 
grazing  angles,  a  modified  two  component  model  was  presented  by  Rosenberg  [2014]  to  model 
the  Ingara  data.  This  led  to  a  simple  model  for  temporal  decorrelation  as  a  function  of  the  wind 
speed  and  wave  height.  This  final  empirical  model  is  presented  in  Section  2.4  with  details  on  how 
to  account  for  platform  motion  and  antenna  scanning. 

2.1  Trials  background 

The  trial  data  was  obtained  with  Ingara  on  two  separate  occasions  and  at  two  distinctly  different 
regions.  The  first  ‘sea-clutter’  trial  was  conducted  in  2004  (SCT04)  in  the  Southern  Ocean 
approximately  100  km  south  of  Port  Lincoln,  South  Australia  [Crisp,  Stacy  &  Goh  2006].  The 
site  chosen  was  at  the  edge  of  the  Australian  continental  shelf  where  there  was  little  chance  of 
shallow  water  affecting  the  wave  field.  During  the  trial,  ocean  backscatter  was  collected  for  a 
range  of  different  geometries  on  eight  separate  days  with  different  ocean  conditions.  The  second 
‘maritime  surveillance’  trial  was  conducted  in  2006  (MAST06)  in  littoral  and  open  sea 
environments  near  Darwin  in  the  Northern  Territory.  Again,  data  was  collected  for  a  range  of 
different  geometries  and  ocean  conditions.  In  this  trial,  a  total  of  four  days  data  were  collected: 
two  of  the  days  were  in  the  littoral  zone  approximately  25  km  north  of  Darwin  and  the  other  two 
days  were  in  the  open  ocean  approximately  200  km  west  of  Darwin. 

During  the  ocean  backscatter  collections,  Ingara  was  operated  in  the  circular  spotlight-mode. 
Figure  1  shows  the  collection  geometry  for  this  mode,  where  the  aircraft  flies  a  circular  orbit  in 
an  anti-clockwise  direction  (as  seen  from  above)  around  a  nominated  point  of  interest,  while  the 
radar  beam  is  continuously  directed  toward  this  point.  Radar  echo  data  is  continuously  collected 
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during  the  full  360°  orbit,  with  the  instantaneous  pulse  repetition  frequency  (PRF)  appropriately 
adjusted  to  maintain  a  constant  spatial  pulse  separation  between  pulse  transmission  positions. 
Once  collected,  the  echo  data  may  be  processed  either  immediately  (in  real-time)  or  subsequently 
(off-line)  to  produce  either  range-compressed  profiles  or  spotlight  synthetic  aperture  radar  (SAR) 
images  of  the  scene  at  various  azimuth  angles.  Further,  owing  to  the  continuous  nature  of  the  data 
collection,  the  images  can  be  formed  at  any  desired  azimuth  look  direction.  Each  collection  of 
data  in  this  mode  is  referred  to  as  a  ‘run’  and  there  may  be  several  complete  orbits  in  a  single  run. 


Figure  1:  Circular  spotlight  mode  collection  geometry  used  for  both  Ingara  sea  clutter  trials. 

In  order  to  examine  the  effect  of  grazing  angle  on  ocean  backscatter,  runs  were  made  with  different 
altitude  and  orbit  diameters.  For  both  the  SCT04  and  MAST06  trials,  data  was  collected  at  the 
centre  of  the  spotlight  for  the  nominal  grazing  angles  of  15°  to  45°  in  5°  increments.  Owing  to 
the  finite  beam  width  of  the  radar,  its  footprint  on  the  ocean  surface  has  a  significant  range  extent. 
This  means  that  the  grazing  angle  varies  across  the  footprint.  It  follows  that,  with  appropriate 
range  compression  and  data  processing,  the  variation  in  backscatter  with  grazing  angle  across  the 
range  extent  of  the  radar  beam  footprint  can  be  measured.  In  this  way,  backscatter  measurements 
for  most  grazing  angles  in  the  range  from  15°  to  45°  could  be  extracted  from  the  data. 

The  aircraft  speed  was  approximately  200  knots  and  so  a  1.5  nautical  mile  (NM)  radius  orbit  took 
approximately  3  minutes  while  a  1.9  NM  radius  orbit  took  about  4  minutes.  The  total  collection 
across  all  grazing  angles  took  approximately  90  minutes.  It  is  reasonable  to  assume  that  over 
such  short  time  intervals,  the  ocean  surface  conditions  are  relatively  unchanged  and  that  mean 
backscatter  variations  are  mostly  due  to  the  changing  imaging  geometry  rather  than  changing 
ocean  conditions.  Nevertheless,  it  is  possible  that  wind  gusts  and  changes  of  wind  strength  and 
direction  may  have  affected  the  measurements.  Table  1  shows  the  wind  and  wave  ground  truth  for 
the  data  used  in  this  report. 

Finally,  using  the  geometry  in  Figure  1 ,  the  azimuth  resolution  can  be  calculated  approximately 
by  R4>3dB,  where  R  is  the  slant  range  and  (fidB  is  the  measured  two-way  azimuth  antenna  3  dB 
beam  width  equal  to  1.02°  for  the  horizontal  transmit,  horizontal  receive  (HH)  channel  and  0.99° 
for  the  vertical  transmit,  vertical  receive  (VV).  Table  2  shows  the  geometry  for  collections  at  a 
range  of  grazing  angles  and  slant-ranges  for  a  beam  width  of  1°.  Over  these  nominal  parameters, 
the  average  azimuth  resolution  is  62.7  m.  Flowever  since  each  collection  spans  a  range  of  grazing 
angles,  the  actual  azimuth  resolution  will  always  differ  slightly. 
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Table  1:  Wind  and  wave  ground  truth.  Directions  are  “from”  not  “to”. 


Trial 

Flight 

Date 

Wind 

Wave 

Speed 

(m/s) 

Direction 

(°) 

Height 

(m) 

Direction 

(°) 

Period 

(s) 

SCT04 

F33 

9/8/04 

10.2 

248 

4.9 

220 

12.3 

SCT04 

F34 

10/8/04 

7.9 

248 

3.5 

205 

11.8 

SCT04 

F35 

11/8/04 

10.3 

315 

2.6 

210 

10.4 

SCT04 

F36 

12/8/04 

13.6 

0 

3.2 

293 

OO 

OO 

SCT04 

F37 

16/8/04 

9.3 

68 

2.5 

169 

9.7 

SCT04 

F39 

20/8/04 

9.5 

315 

3.0 

234 

11.4 

SCT04 

F40 

24/8/04 

13.2 

22 

3.8 

254 

12.2 

SCT04 

F42 

27/8/04 

8.5 

0 

4.3 

243 

12.5 

MAST06 

F2 

17/5/06 

8.5 

115 

0.62 

112 

3.1 

MAST06 

F4 

19/5/06 

3.6 

66 

0.25 

35 

2.6 

MAST06 

F8 

23/5/06 

3.5 

83 

0.41 

46 

4.0 

MAST06 

F9 

24/5/06 

10.2 

124 

1.21 

128 

4.6 

Table  2:  Nominal  geometric  parameters  of  circular  spotlight-mode  collections. 


Grazing  angle  (°) 

Altitude  (m) 

Radius  (NM  /  m) 

Slant  range  (m) 

Azimuth  res.  (m) 

15 

932 

1.9/3519 

3643 

63.54 

20 

1353 

1.8/3334 

3548 

62.80 

25 

1522 

1.8/3334 

3679 

63.97 

30 

1711 

1.6/2963 

3421 

59.72 

35 

2073 

1.6/2963 

3617 

63.11 

40 

2314 

1.5/2778 

3626 

63.10 

2.2  Clutter  and  noise  mean  power 


Simulating  radar  performance  requires  accurate  models  for  the  mean  power  of  the  clutter,  noise 
and  target.  To  measure  the  performance,  the  target  power  is  typically  varied  with  respect  to  the 
interference  or  clutter  plus  noise  power.  The  mean  clutter  power,  pc  is  determined  by  the  radar 
range  equation: 


PtG2  \2ca  CTPB 

(4vr)3JR4LaLs 


(1) 


where  Pt  is  the  transmit  power,  G  is  the  one-way  gain  on  transmit  and  receive,  Ac  is  the  radar 
wavelength,  R  is  the  slant  range,  ac  is  the  sea-clutter  RCS,  La  and  Ls  represent  the  atmospheric 
and  system  losses  respectively  and  pulse  compression  adds  a  gain  given  by  the  pulse  length  - 
bandwidth  product,  TpB.  The  thermal  noise  power  in  the  radar  can  be  defined  as, 


Pn  =  kboltzT0BFn  (2) 

where  kbobz  is  Boltzmann’s  constant,  To  is  the  reference  temperature  and  Fn  is  the  radar  noise 
figure.  The  next  section  gives  further  detail  of  the  parameters  in  Equation  1.  Section  2.2.2  then 
presents  the  Imaging  Radar  Systems  Group  (IRSG)  mean  backscatter  model  which  captures  the 
mean  backscatter  as  a  function  of  azimuth,  grazing  and  wind  speed.  These  are  then  all  brought 
together  with  a  model  for  the  CNR  in  Section  2.2.3. 
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2.2.1  Parameter  description 

There  are  many  parameters  in  Equation  1  which  are  required  to  accurate  represent  the  CNR.  The 
values  chosen  for  this  analysis  arc  now  described  in  more  detail: 


Transmit  power,  Pt  The  peak  power  of  a  typical  airborne  radar  is  20  kW. 

Radar  wavelength,  Ac  For  an  X-band  radar  with  a  centre  frequency,  fc  of  10  GHz,  the 
wavelength  is  3  cm. 

Antenna  gain,  G  The  one-way  antenna  gain  is  related  to  the  wavelength  and  the  effective  area  of 
the  antenna  aperture.  For  the  modelling  in  this  report,  an  antenna  gain  of  35  dB  has  been 
used. 


Slant  range,  R  The  slant  range  is  related  to  the  altitude  of  the  radar,  h,  the  grazing  angle,  0  and 
the  mean  radius  of  the  earth,  Re  =  6.4  x  103  km.  The  earth  radius  is  typically  multiplied 
by  4/3  to  account  for  the  average  effect  of  atmospheric  refraction,  giving  an  effective  radius 
E  =  4/3 Re.  The  slant  range  to  the  sea  surface  is  obtained  by  solving  equation  (12.10)  of 
[Ward,  Tough  &  Watts  2013]  for  R: 

R  =  \JE 2  sin2  9  +  2 Eh  +  K2  -  E  sin  9 
^  h(l  +  h/E) 
sin  9  +  h/E 

The  flat  earth  approximation  is  R  «  hj  sin  9,  which  is  quite  accurate  unless  the  grazing 
angle  is  small.  For  example,  if  the  radar  altitude  is  10  km  (33.1  thousand  feet)  and  the 
grazing  angle  is  30°,  the  flat  earth  approximation  for  the  slant  range  is  20  km,  whereas  the 
more  accurate  calculation  gives  a  result  just  35  m  less  than  this. 

Sea-clutter  RCS,  ac  The  radar  cross  section1  can  be  described  in  terms  of  the  mean  backscatter 
a o  within  a  given  resolution  cell  with  area  Acen, 

ac  =  o-0Aceii.  (4) 


The  mean  backscatter  is  modelled  using  the  IRSG  model  described  below,  while  the  ground 
plane  area,  Ace\\  is  given  by  the  slant  range  and  azimuth  resolutions,  Ax,  Ay  projected  onto 
the  ground, 

Aceii  —  AxAy  sec  9  (5) 

with  range  and  azimuth  resolutions  given  by, 


A*  =  — , 

*  2  B' 

Ay  «  R4>3dB 


(6) 


where  vc  is  the  speed  of  light.  For  a  bandwidth  of  200  MHz,  the  range  resolution  is  0.75  m. 
For  the  azimuth  resolution,  with  a  beam  width  of  1°,  the  resolution  varies  according  to  the 
grazing  angle.  For  the  Ingara  flight  trials,  the  azimuth  resolution  was  approximately  63  m 
(see  Table  2). 

'Note  that  the  propagation  factor  is  inherently  included  in  this  definition.  Atmospheric  losses  have  been  included  in 
the  loss  term,  La- 
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Atmospheric  loss,  La  The  atmospheric  loss,  La  describes  all  the  losses  in  the  atmosphere.  A  loss 
of  2  dB  has  been  chosen  for  this  analysis. 

System  loss,  Ls  The  system  loss,  Ls  is  a  bulk  term  which  describes  all  the  losses  in  the  radar 
system.  A  loss  of  5  dB  has  been  chosen  for  this  analysis. 

Pulse  compression  gain,  TpB  Pulse  compression  adds  a  gain  to  the  radar  range  equation  based 
on  the  pulse  length,  Tp  and  the  system  bandwidth.  For  a  pulse  width  of  5//.s  and  a  bandwidth 
of  200  MHz,  the  gain  is  1000. 

Temperature,  To  The  ambient  temperature  is  set  to  290  K. 

Noise  figure,  Fn  The  system  noise  figure  is  a  measure  of  the  degradation  of  the  signal  to  noise 
ratio  caused  by  components  in  the  RF  signal  chain.  A  typical  value  is  4  dB. 

In  addition  to  these  parameters,  the  models  for  the  mean  backscatter  and  Pareto  shape  require  the 
wind  speed  and  significant  wave  height  as  inputs.  For  modelling  puiposes,  it  is  desirable  to  specify 
a  single  parameter  to  represent  these  sea  conditions.  The  Douglas  sea-state,  S  can  be  related  to 
both  wind  speed  and  significant  wave  height  using  the  following  relationships  [Technology  Service 
Corporation  1990]: 


U  =  3.2  S0-8, 
H1/3  =  0.024 U2. 


(V) 

(B) 


2.2.2  IRSG  mean  backscatter  model 

There  are  a  number  of  existing  models  for  the  mean  backscatter,  although  only  a  few  cover  the 
medium  grazing  angle  region.  The  two  most  suited  are  the  Technology  Services  Corporation 
(TSC)  and  Masuko  models,  [Technology  Service  Corporation  1990,  Spaulding,  Horton  &  Pham 
2005,  Masuko  et  al.  1986],  however  they  are  both  unable  to  capture  all  the  variations  with  the 
weather  conditions  and  azimuth  and  grazing  angles.  This  led  to  a  new  empirical  model  presented 
by  [Crisp  et  al.  2008]  known  as  the  imaging  radar  systems  group  (IRSG)  model.  The  model  is 
applied  independently  to  each  polarisation  in  the  log  domain  and  uses  relationships  taken  from 
[Ulaby,  Moore  &  Fung  1982,  Shimada,  Kawamura  &  Shimada  2003].  It  includes  a  power  law 
model  for  variations  in  the  wind  speed,  a  linear  model  for  the  grazing  and  a  two  element  Fourier 
series  to  model  azimuth.  The  principle  behind  this  model  is  that  the  small  scale  waves  which 
are  responsible  for  Bragg  scattering  from  the  ocean  surface  are  directly  dependent  on  the  surface 
wind  speed.  The  stronger  the  wind  speed,  the  larger  the  capillary  waves  and  hence  the  larger  the 
backscatter.  The  model  is  valid  over  20°-45°  in  grazing,  360°  in  azimuth  and  wind  speeds  from 
3-14  m/s. 

The  first  component  of  the  IRSG  model  captures  the  variation  in  backscatter  with  grazing  angle 
0  and  wind  speed  U  and  is  defined  for  the  upwind,  downwind  and  crosswind  directions.  If  a'0 
represents  the  mean  backscatter  in  one  of  these  azimuth  directions,  then  the  relationship  is  given 
by 


a'0(6,  U)  =  a0  +  aid  +  a2  log10(U). 


(9) 
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where  the  three  sets  of  coefficients  ao,  «i  and  02  have  been  determined  by  model  fitting  over  the 
12  days  of  the  Ingara  data.  The  second  stage  is  to  combine  these  three  backscatter  results  into  new 
coefficients, 


b0(e,  U )  =  [o>u{6,  U)  +  cr'd(9,  U)  +  2 </c(0,  U)]  /4, 

W  U)  =  [(o'Jfl,  U)  -  </d(0,  U)\  /2,  (10) 

62(0,U)  =  6o(0,U)-a'(0,C/), 

where  cr^(-),  crd(-)  and  ay  ( • )  are  the  upwind,  downwind  and  cross  wind  values  of  a'0(-).  The  final 
step  takes  these  new  coefficients  and  includes  the  azimuth  variation,  </>,  in  a  two  element  Fourier 
series, 

<r0(<£,  0,  U )  =  b0(e ,  U)  +  61(0,  U)  cos  <j>  +  62(6>,  U)  cos  (2<£)  (1 1) 

Implementing  the  model  requires  calculating  Equation  9  for  each  wind  direction  using  three  sets  of 
a  coefficients  for  each  azimuth  direction.  The  mean  backscatter  is  then  determined  by  substituting 
the  angle  specific  backscatter  values  into  Equation  10  and  then  using  that  result  in  Equation  11. 
The  9  coefficients  required  to  implement  the  model  are  given  in  Table  3  for  each  polarisation. 
Figure  2  shows  an  example  of  the  F35  data  set  with  the  wind  speed  taken  from  Table  1.  Clearly 
the  trends  have  been  well  captured  in  the  model  with  a  sinusoidal  variation  across  azimuth  and 
increasing  mean  backscatter  as  the  grazing  angle  increases. 

Table  3:  Mean  backscatter  coefficients. 


Azimuth 

ao 

at 

a2 

Upwind 

-60.03 

23.39 

22.65 

HH 

Downwind 

-67.80 

28.58 

23.92 

Crosswind 

-67.09 

23.12 

24.71 

Upwind 

-66.65 

9.72 

25.46 

HV 

Downwind 

-68.74 

14.17 

24.47 

Crosswind 

-73.09 

14.85 

26.66 

Upwind 

-50.18 

12.41 

25.15 

VV 

Downwind 

-50.16 

12.30 

23.92 

Crosswind 

-52.60 

12.30 

22.09 

2.2.3  Clutter  to  noise  ratio 

To  demonstrate  the  modelling  of  these  parameters,  the  CNR  has  been  calculated  using  parameters 
representative  of  the  F35  flight  with  the  mean  backscatter  determined  from  the  IRSG  model. 
Figure  3  shows  the  result  which  demonstrates  a  good  match  with  appropriate  variation  in  both 
azimuth  and  grazing. 


2.3  Probability  distributions 

Statistical  models  of  the  PDF  are  used  to  model  the  spread  of  received  backscatter  values.  A 
number  of  them  have  been  used  to  model  sea-clutter  in  the  past  with  varying  degrees  of  success 
[Dong  2006,  Rosenberg,  Crisp  &  Stacy  2010fi,  Weinberg  201  In,  Rosenberg  &  Bocquet  2015, 
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Mean  backscatter 


Azimuth  (deg) 
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Mean  backscatter  comparison  with  varying  azimuth  and  grazing  using  the  IRSG  model, 
left  -  Ingara  F35  data,  right  -  model  output  with  wind  speed,  U  =  10.3  m/s. 
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Figure  3:  CNR  comparison  with  varying  azimuth  and  grazing:  left  -  Ingara  F35  data,  right  - 
model  using  parameters  matched  to  data  and  mean  backscatter  determined  from  the 
IRSG  model. 
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Rosenberg,  Watts  &  Bocquet  2014].  The  first  model  considered  is  the  K-distribution  which  is 
used  to  characterise  the  bulk  of  the  PDF.  To  understand  the  development  of  the  K-distribution, 
consider  a  radar  receiving  in-phase  and  quadrature  data  from  an  external  clutter  source  with  its 
amplitude  defined  by  Gaussian  statistics  with  zero  mean  and  variance,  x.  In  addition,  thermal 
noise  from  within  the  radar  will  add  a  component  pn  which  is  included  by  offsetting  the  variance 
x.  The  PDF  of  the  clutter  and  noise  is  then  given  by, 


PE(E\x,pn ) 


1 

.  - exp 

a/27 T(x+pn) 


E 2 

2(x+pn) 


(12) 


with  the  envelope  of  these  components,  y 


,  then  determined  by  a  Rayleigh  PDF, 


Py(y\x,Pn) 


y 


X+Pn 


exp 


2(x  +Pn). 


(13) 


In  the  target  detection  analysis,  the  received  pulses  are  often  converted  to  power  (square  law), 
z  =  y1  and  then  summed  to  improve  detection  performance.  Consider  the  sum  of  M  exponential 
random  variables, 

M 

Z  =Y2zm-  U4) 

?7I=1 

If  the  returns  from  successive  pulses  are  uncorrelated,  the  received  power  is  then  described  by  a 
gamma  PDF, 

r^M—1 

Pl'-(z'I)=(.>+p.)"r(M)elp 


X+Pn 


(15) 


2.3.1  Compound  distributions 

The  distribution  in  Equation  15  models  the  speckle  component  of  the  sea-clutter.  To  include  the 
texture  component  which  modulates  the  speckle,  a  compound  distribution  can  be  used, 

noo 

P{Z)  =  /  Pz]x(Z\x)Px(x)dx  (16) 

J  o 

where  Px(x)  is  the  distribution  of  the  texture  component.  To  achieve  a  K-distribution,  the  texture 
is  also  gamma  distributed 


Px(x)  =  exP  i~bx\  >  U  b  >  0 


(17) 


with  shape  given  by  v  and  scale,  b  =  v /pc  related  to  the  mean  power,  pc.  For  the  Pareto 
distribution,  the  texture  has  an  inverse  gamma  distribution 


Px{x)  =  f^x  a  1  exp  [  cf xj  ,  a  >  1,  c  >  0 
1  (aj 


(18) 


where  a  is  the  shape  and  c  =  pc(a  —  1)  is  the  scale.  In  general,  numerical  integration  is  required  to 
evaluate  the  compound  distributions.  However,  if  thermal  noise  is  not  present  (pn  =  0),  analytic 
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forms  of  the  distributions  can  be  derived.  With  a  gamma  distributed  texture,  the  K-distribution  is 
obtained, 


2  M+v  1  /  , - \ 

P(Z)  =  -  ( bZ )—  — — - — —Ku_m  (  2 VbZ)  (19) 

V  ;  Z  V  ’  r(M)I»  V  ) 

where  is  the  modified  Bessel  function  of  the  second  kind  with  order  v—M.  If  the  texture 

is  inverse  gamma  distributed,  the  general  solution  is  a  compound  gamma  distribution 


pm=  ZM~1caT(M  +  a) 

(c  +  z)M+“T(M)r(ay 


a  >  1,  c  >  0 


(20) 


which  reduces  to  the  generalised  Pareto  type  II  or  Lomax  distribution  for  M  =  1, 


P(Z) 


aca 

( c+Z)l+a ' 


(21) 


The  result  in  Figure  4  shows  the  time  domain  PDF  and  the  complementary  cumulative  distribution 
function  (CCDF)  for  the  F35  dataset  (in  intensity)  with  the  Pareto  distribution  fit  overlayed  in 
black.  For  comparison,  the  equivalent  fit  for  the  K  distribution  is  shown  in  red  and  is  clearly 
under-fitting  the  data,  particularly  in  the  F1H  channel.  Thermal  noise  is  assumed  for  both  models 
with  the  recently  published  z  log  z  shape  estimation  method  used  to  estimate  the  shape  values 
[Bocquet  2015]. 


PDF  CCDF 


Relative  backscatter  (dB) 


Relative  threshold  (dB) 


Figure  4:  Time  domain  intensity  distributions  -  upwind  30°  grazing.  Data  in  blue,  K-distribution 
fit  in  red,  Pareto  distribution  fit  in  black. 
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2.3.2  Shape  model 

Until  recently,  the  only  model  suitable  for  modelling  the  K-distribution  shape  was  for  low  grazing 
angles  [Ward,  Tough  &  Watts  2013].  Using  the  Ingara  data  set,  there  have  now  been  three  iterations 
of  models  designed  to  model  the  shape  at  medium  grazing  angles  with  all  three  applied  separately 
to  each  polarisation.  The  first  was  designed  to  model  only  the  geometry  for  an  individual  day 
[Rosenberg,  Crisp  &  Stacy  2009],  This  was  then  extended  by  Crisp,  Rosenberg  &  Stacy  [2015] 
to  include  extra  components  for  the  environmental  conditions  and  resolution.  These  are  combined 
into  a  global  model  comprising  24  coefficients  at  the  finest  resolution  (0.75  m  range  x  63  m 
azimuth).  The  third  iteration  is  a  variation  designed  to  better  model  other  parameters  of  interest 
which  have  a  stronger  component  due  to  the  swell  [Rosenberg  et  al.  2015].  It  was  also  applied 
to  model  the  Pareto  shape  and  offers  a  slightly  better  fit  to  the  estimated  K-distribution  shape. 
Consequently,  this  third  iteration  is  used  in  the  report.  The  validity  of  the  model  parameter  bounds 
are  determined  by  the  input  data.  This  includes  a  CPI  of  64  pulses,  a  resolution  of  0.75  m  range  x 
63  m  azimuth,  all  azimuth  angles  and  linear  polarisations,  sea-states  2-6  and  grazing  angles  from 
15°  -  45°. 

If  X  is  a  generic  parameter  to  be  modelled,  then  the  geometric  relationship  can  be  written  as  the 
weighted  sum  of  two  components  which  represent  the  wind  and  swell  directions.  Assuming  the 
data  has  been  rotated  so  upwind  is  now  at  0°, 

X(0,  4>)  =  a'o07  [l  +  a\  cos  <f>  +  a2  cos(2 (f>)  +  a'3  cos (</>  —  ip)  +  a'4  cos(2(i ft  —  ^>))]  (22) 

where  ^  is  the  swell  direction  relative  to  upwind  and  7,  a'0, ...  ,a'4  are  the  model  coefficients. 
Representing  the  ocean  conditions  is  difficult  and  the  shape  model  uses  only  the  wind  speed  and 
the  significant  wave  height  -ffi/3, 

Y  =  b,0  +  b'1log10(U)  +  b,2H1/3  (23) 

where  b'Q,  b\ ,  b2  are  the  model  coefficients.  To  relate  these  two  models,  the  coefficients  in  Equation 
22  must  be  altered  so  they  are  independent  of  grazing  angle.  This  is  achieved  by  introducing  a 
normalisation  factor,  6q  and  then  redefining  Equation  22  as 

/  Q  y 

X(8,<f>)=[  —  \  [ao  +  Oi\  cos  4>  +  «2  cos(2cj))  +  CK3  cos(0  —  ip)  +  04  cos(2(c/)  —  y)))]  (24) 

where  the  new  coefficients  are  related  by 

ao  =  a'o0Q,  ai  =  a'oa^Q,  ...,  07  =  a'ga^Q.  (25) 

Implementing  the  model  involves  matching  7  and  each  of  the  a  coefficients  to  Equation  23.  This 
in  turn  requires  6  sets  of  the  b'  parameters  or  18  coefficients  in  total  for  each  polarisation.  These 
model  coefficients  for  K  and  Pareto  distributions  are  presented  in  Table  4.  Once  the  new 
coefficients  in  Equation  25  are  evaluated,  they  can  be  substituted  into  Equation  24  with  the 
normalising  grazing  angle  6q  set  to  30°  to  derive  the  final  model  value. 

Finally,  it  has  been  observed  that  shape  values  less  than  2  can  result  in  bad  fits  to  the  Ingara  data. 
To  enforce  this  lower  bound  with  the  new  shape  model,  the  model  coefficients  have  been  derived 
for  a  modified  shape,  a  —  2,  with  the  final  shape  value  then  increased  by  2  to  compensate.  Figure 
5  shows  the  variation  of  the  Pareto  shape  over  grazing  and  azimuth  for  the  F35  run  with  a  model 
created  with  the  same  input  parameters.  Clearly,  the  model  captures  the  main  trends  quite  well. 
These  include  a  include  a  linear  increase  with  grazing  and  sinusoidal  variation  in  azimuth  with  a 
peak  in  the  upwind  and  downwind  directions. 
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Figure  5:  Pareto  shape  parameter  model  comparison  with  varying  azimuth  and  grazing:  left  - 
shape  estimated  from  F35  data  set,  right  -  model  output  with  wind  speed,  U  =  10.3  m/s, 
significant  wave  height,  II  \  /^  =  2.6  m,  wind  swell  angle,  =  0°. 


2.4  Temporal  correlation 


To  accurately  capture  the  correlation  due  to  the  speckle,  platform  motion  and  antenna  scanning,  a 
model  is  required  which  relates  both  the  sea  conditions  and  the  relevant  radar  parameters.  In  the 
following  sections,  different  Gaussian  models  are  used  to  model  each  of  the  components  in  the 
frequency  domain  with  their  overall  spread  spectrum  given  as  their  convolution: 


C(f)  = 


1 


y/2- 7T 


a 


exp 


(/  -  fo/ 

2cr2 


where  /o  is  the  mean  and  the  total  variance  is  the  sum  of  the  component  variances, 


a'2  ~  as  +  apm  +  ’-’scan • 


(26) 


(27) 


The  combined  temporal  auto-correlation  function  (ACF)  is  then  found  by  taking  the  Fourier 
transform 


C(t)  =  exp(— j27r/0f) exp  ( 


with  the  overall  decorrelation  time  given  by 

T'  = 


y/2- 


(28) 


(29) 


TUT 
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Table  4:  K  and  Pareto  distribution  shape  coefficients. 


K  shape 

Pareto  shape 

b'o 

b'2 

b'o 

b\ 

b'2 

7 

11.52 

-8.72 

-0.11 

13.4 

-10.15 

-0.064 

-1.37 

4.18 

-0.078 

-0.93 

3.00 

-0.11 

HH 

a\ 

-0.38 

0.11 

0.050 

-0.39 

0.15 

0.065 

OL2 

-0.14 

0.80 

0.024 

-0.31 

0.89 

-0.028 

«3 

-0.018 

-0.19 

0.043 

0.32 

-0.49 

0.050 

-0.51 

0.68 

0.019 

-0.62 

1.01 

0.00037 

7 

4.83 

-3.33 

0.054 

8.34 

-6.64 

0.061 

«o 

-5.78 

11.52 

-0.39 

-5.68 

10.75 

-0.40 

HV 

a\ 

-4.17 

6.38 

-0.24 

-3.55 

5.36 

-0.21 

OL2 

-0.89 

3.32 

-0.033 

-0.52 

2.72 

-0.21 

OL  3 

1.98 

-2.59 

-0.030 

2.96 

-4.35 

0.080 

a  4 

-0.16 

0.12 

0.40 

0.41 

0.48 

0.22 

7 

0.28 

-0.38 

0.11 

0.28 

-0.39 

0.11 

ao 

14.06 

19.11 

-4.06 

13.43 

18.78 

-4.06 

VV 

a\ 

-6.00 

3.80 

0.60 

-7.37 

5.61 

0.64 

OL2 

0.93 

8.25 

-0.80 

1.41 

7.29 

-1.03 

«3 

-0.72 

0.15 

0.16 

2.47 

-2.49 

0.095 

a  4 

-3.91 

5.56 

0.40 

-10.87 

16.67 

-0.41 

2.4.1  Clutter  width 


As  the  sea  state  becomes  greater,  the  wind  speed,  wave  height,  wave  length  and  the  Doppler 
bandwidth  all  increase.  Regarding  the  latter,  Hicks  et  al.  [1960]  have  shown  experimentally  that 
the  Doppler  spread  at  high  grazing  angles  is  proportional  to  the  wave  height  divided  by  the  wave 
period.  Previous  analysis  of  the  bi-modal  Doppler  spectrum  [Rosenberg  2014]  found  only  minor 
trends  in  azimuth  and  grazing  angle  and  hence  the  model  is  only  a  function  of  wind  speed  and 
wave  height.  The  model  represents  the  temporal  decorrelation  which  is  measured  at  the  point 
where  the  absolute  value  of  the  ACF  decays  to  1/e.  It  is  described  by  three  sets  of  coefficients, 
one  for  each  polarisation.  Note  that  this  model  only  considers  the  real  part  of  the  correlation  and 
hence  assumes  that  the  Doppler  spectrum  is  centred  at  0  Hz.  The  model  is  given  by 

T  =  (po  +PiU  +  p2H1/3)  x  1(T3  (30) 


with  the  model  coefficients  given  in  Table  5.  Figure  6  shows  the  model  fit  for  the  two  component 
model  where  the  trend  is  that  the  mean  decorrelation  times  arc  longer  when  the  wind  speed  and 
wave  height  are  smaller.  The  RMS  errors  for  the  HH  and  VV  channels  are  1.1  and  0.74  ms 
respectively.  The  HV  decorrelation  time  is  modelled  as  the  average  of  the  HH  and  VV  model 
outputs. 

A  Gaussian  model  can  then  be  used  to  represent  the  correlation, 


A(t)  =  exp 


(31) 
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Table  5:  Decorrelation  model  parameters 


HH 

VV 

Po 

16.67 

17.13 

Pi 

-0.68 

-0.40 

P2 

-0.87 

-0.91 

Figure  6:  Mean  temporal  decorrelation  as  a  function  of  wind  speed  and  wave  height.  Lefi  -  FIFl, 
Right  -  W. 

with  the  spectrum  found  by  taking  a  Fourier  transform 

A(f)  =  s/nT  exp  (— 7r2/2T2) 

=  ^exp(^)  (32) 

where  <rs  =  1/(V2ttT). 


2.4.2  Platform  motion 


Antenna  beampatterns  are  typically  defined  in  an  angular  domain  with  specifications  for  the 
minimum  sidelobe  level  and  two-way  3  dB  azimuth  beam  width.  To  simplify  the  modelling,  we 
will  assume  a  Gaussian  beam  pattern  with  azimuth  angle,  cj) o, 


G(cj>  0)  =  exp 


(21n2)03dB 


(33) 


Assuming  that  the  elevation  beam  pattern  does  not  change  significantly  over  the  clutter  patch,  the 
Doppler  frequency  can  be  related  to  the  beam  pattern  by 


/ 


2vpfa  |  sin  cj) b  |  cos  9e i 

A 


(34) 


where  fa  is  the  antenna  azimuth  direction  relative  to  the  platform  velocity  and  9e\  is  the  elevation 
angle.  The  azimuth  angle  can  then  be  related  to  the  Doppler  frequency  and  hence  the  platform 
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velocity  by  substituting 


into  Equation  33,  giving 


where 


2vp\  sin^>b|  cos  9e i 


G(/)R<exp(^) 

_  1  vP\  sin  <jib|03dB  cos  (9ei 

Pm  “  ^2!n2  A 

_  n  a-vp\  sin  0b|(/>3dB  cos  9e\ 


(35) 


(36) 


(37) 


2.4.3  Antenna  scanning 

Although  the  Ingara  radar  did  not  scan  during  the  sea-clutter  collections,  antenna  scanning  can  be 
modelled  using  the  relationship  in  Skolnik  [2008], 

\/2  In  2  (/>scan 

^scan  =  A  7 

<P3dB 

»  0.19^^  (38) 

<P3dB 

where  </>scan  is  the  antenna  rotation  rate  in  rad/s. 
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3  Detection  probability  framework 

The  problem  of  target  detection  is  commonly  presented  with  two  hypotheses,  Hq  -  clutter  plus 
noise  is  present  in  the  radar  return  and  H\  -  a  target  is  present  with  the  clutter  and  noise.  Figure 
7  shows  the  case  for  non-coherent  detection  where  two  histograms  are  shown  which  correspond 
to  these  hypotheses.  The  premise  of  target  detection  is  to  select  an  acceptable  probability  of  false 
alarm,  Pfa  and  determine  the  threshold,  Y  based  on  the  clutter  plus  noise  model.  Any  backscatter 
response  above  the  threshold  is  then  declared  a  potential  target.  If  a  detector  can  measure  or 
estimate  the  clutter  distribution  parameters  and  set  the  corresponding  threshold,  then  the  detector 
can  be  considered  CFAR.  The  Pfa  and  the  probability  of  detection,  If  can  be  defined  as 


P&(Y)  =  Prob (Z  >  Y\H0)  = 

/•oo 

J  Pc(Z)dZ 

(39) 

Pd(Y)  =  Prob(Z  >  Y\Hi)  = 

/»oo 

J  Pt(Z)dZ 

(40) 

where  PC(Z )  represents  the  clutter  plus  noise  PDF  and  If ( Z )  the  target  plus  clutter  and  noise  PDF. 
Also  in  this  figure,  Pm  denotes  the  miss  probability. 


Figure  7:  Detection  probabilities  showing  distributions  for  the  clutter  plus  noise  and  the  target 
plus  clutter  and  noise.  The  threshold  determines  the  three  regions:  P,„,  Pfa  and  If. 

The  theory  of  optimal  target  detection  is  well  summarised  in  Ward,  Tough  &  Watts  [2013,  Chapter 
10].  The  principle  of  finding  an  optimal  Ff  for  a  fixed  f’a  is  known  as  the  Neyman-Pe  arson 
criterion  and  can  be  written  as  the  ratio  of  the  interference  and  target  plus  interference  PDFs 

k^=W)>Y-  <41) 

Expressions  for  compound  PDFs  are  generally  quite  complex,  particularly  when  correlation  is 
included  and  the  signals  are  treated  coherently.  Also  CFAR  control  remains  difficult  when 
estimates  of  the  covariance  are  required  and  the  sea-clutter  is  non- stationary.  This  has  led  to 
non-coherent  schemes  which  are  often  sub-optimal  and  may  or  may  not  be  CFAR.  These 
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practical  detectors  are  typically  easier  to  compute  and  require  adjustment  of  the  threshold  to 
maintain  a  constant  false  alarm  rate.  Some  of  the  latter  include  the  cell- averaging,  single-sided 
and  ordered  statistic  CFAR  [Ward,  Tough  &  Watts  2013,  Chapter  13].  This  report  only  considers 
detection  with  a  constant  or  fixed  threshold  and  with  non-coherent  integration  as  the  mathematics 
are  more  tractable.  We  assume  that  a  threshold  is  applied  after  pulse-to-pulse  integration  and 
there  is  a  CFAR  system  to  adjust  the  threshold  for  a  specified  Pfa.  The  exact  analytic  expressions 
then  represent  the  performance  results  which  would  be  expected  if  the  clutter  parameters  are 
known  and  are  useful  as  a  benchmark  for  other  detection  approaches. 

This  section  provides  details  of  the  Pfa  and  P(j  calculations  for  a  target  present  in  clutter  and 
noise  (interference).  The  first  Sections  3. 1-3.2  assume  that  the  clutter  and  target  distributions 
are  uncorrelated,  such  as  would  be  the  case  in  a  sufficiently  frequency  agile  radar.  Section  3.1 
considers  a  constant  target  present  in  Gaussian  interference  covering  a  ‘local  area’.  The  solution 
for  the  Pd  is  known  as  a  Marcum  model  and  was  first  described  by  Marcum  [I960].  Flowever, 
the  RCS  of  a  target  will  typically  fluctuate  and  Section  3.2  describes  a  method  for  representing 
this  fluctuation  in  terms  of  generalised  Weinstock  fluctuations,  [Weinstock  1964]  and  Swerling 
models,  [Swerling  I960]. 

Some  level  of  temporal  correlation  will  typically  be  present  on  both  the  target  and  the  speckle 
component  of  the  clutter.  For  the  former,  it  is  commonly  assumed  to  be  either  totally  correlated 
or  uncorrelated  and  fluctuate  around  a  mean  level  with  an  exponential  distribution  (Swerling 
cases  1  and  2)  or  chi-squared  distribution  (cases  3  and  4).  Flowever,  an  early  study  by  Edrington 
[1965]  found  that  radar  echoes  from  a  number  of  aircraft  were  in  the  fact  partially  correlated  and 
that  performance  prediction  models  needed  to  be  modified.  Barton  [1969]  looked  at  the  target 
correlation  and  calculated  an  effective  number  of  independent  looks  based  on  the  correlation  time 
relative  to  the  total  observation  time.  With  this  approach,  the  Pci  can  then  be  determined  using  the 
traditional  solution  for  totally  correlated  targets.  A  subsequent  study  by  Kanter  [1986]  then 
derived  the  exact  Pd  for  exponential  fluctuating  targets  rather  than  with  an  effective  number  of 
looks.  Weiner  [1988]  then  did  the  same  for  chi-squared  fluctuating  targets. 

Regarding  the  speckle  component,  there  are  a  number  of  approaches  in  the  literature  which  can 
be  used  to  account  for  different  levels  of  correlation.  Ideally,  an  approach  should  be  able  to  model 
the  ACF  through  its  representation  of  the  Doppler  spectrum  of  the  received  signal,  which  is 
related  through  a  Fourier  transform.  The  simplest  method  was  presented  by  Farina  et  al.  [1996] 
and  is  based  on  an  improvement  factor  which  can  be  used  to  modify  the  required  SIR  in  the  Pd 
calculation.  However,  it  assumes  a  Gaussian  Doppler  spectrum  and  hence  cannot  be  used  to 
model  more  complex  spectra  with  components  due  to  breaking  waves  and  discrete  spikes.  A 
second  approach  which  is  only  suitable  for  exponential  fluctuating  targets  was  presented  by  Hou 
&  Morinaga  [1988],  This  result  allows  both  partially  correlated  targets  and  speckle  with  exact 
representations  of  their  ACF’s,  but  suffers  from  numerical  instabilities  and  is  not  generalised  for 
all  Swerling  target  models.  The  final  approach  from  the  literature  by  Ward,  Tough  &  Watts 
[2013]  (WTW)  treats  the  speckle,  thermal  noise  and  target  separately.  By  using  an  estimate  of  the 
‘effective’  number  of  looks  for  the  speckle,  it  can  overcome  the  numerical  problems  of  Hu  and 
Morinaga.  While  this  method  is  useful,  it  becomes  very  computationally  intensive  as  the  CNR 
gets  higher  and  fails  completely  when  there  is  only  clutter  present.  This  led  to  a  modification  by 
Rosenberg  &  Bocquet  [2013ft]  which  overcomes  this  limitation  by  allowing  any  positive  number 
of  looks  for  either  the  speckle  or  target.  Thermal  noise  has  also  been  included  in  this  formulation, 
with  the  caveat  that  it  will  be  correlated  with  the  speckle  and  is  therefore  not  as  accurate  at  low 
CNR.  A  tradeoff  is  therefore  proposed,  where  the  WTW  method  is  used  for  low  CNR  and  the 
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modification  used  for  larger  CNR.  These  algorithms  are  described  in  Section  3.3  with  the  last 
Section  3.4  outlining  how  the  models  can  be  extended  to  a  compound  representation  where  the 
local  area  is  now  augmented  by  a  ‘global’  fluctuation.  This  ‘texture’  is  assumed  to  be 
uncorrelated  and  hence  the  results  presented  here  may  not  be  completely  accurate  for 
scan-to-scan  radars  where  a  longer  time  between  looks  may  result  in  a  variation  of  the  texture. 


3.1  Detection  with  a  constant  target 

The  first  model  represents  a  target  with  a  non-fluctuating  or  constant  RCS  present  in  Gaussian 
interference  and  is  known  as  a  Marcum  model  or  Swerling  case  0.  In  the  intensity  domain,  the 
PDF  for  a  single  pulse  is  given  by  a  Rice  distribution,  [Ward,  Tough  &  Watts  2013], 


P{z\x,pn,A)  = 


1 


exp 


V  X  +  pnJ  \x+pnJ 


(42) 


X+Pn 

where  the  target  power  is  given  by  A2  and  Iq(-)  is  a  modified  Bessel  function  of  the  first  kind, 
order  0.  For  the  following  derivations,  it  is  useful  to  define  the  local  received  signal  and  the  target 
power  as  being  normalised  to  the  interference  power  (i.e.  the  local  SIR), 


1 


M  = 


X+Pn 


S  = 


1 

X+Pn 


M 

^  ]  zmi 

(43) 

m=  1 

M 

(44) 

m=  1 


where  x  and  pn  are  the  speckle  and  thermal  noise  means  and  M  is  the  number  of  looks  or  non¬ 
coherent  averages.  The  PDF  describing  the  interference  is  then  given  by, 


uM~1 

p M  =  iW"" 


(45) 


and  the  Pfa  can  be  determined  by  calculating 

Pfa(j\M)  = 


(46) 


p(+)d^ 

_  r  (M,r) 

r(M) 

where  t  is  the  normalised  threshold  related  to  the  ‘true’  threshold  Y 

Y 

t  =  — -  (47) 

X+Pn 

and  r(M)  =  (M  —  1)!  and  T(M,t)  is  the  upper  incomplete  gamma  function  which  can  be 
evaluated  for  integer  M  using  the  following  series  [Gradshteyn  &  Ryzhik  1994,  §3.351/2  and 
§8.352/2]: 

poo 

M—l  —L 


T(M,t)  =  J  pM~Le-^dpi 
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=  T(M)e~T  £  " 


k= 0 
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(48) 
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For  the  total  received  signal,  the  PDF  becomes  a  multi-look  Rice  distribution, 

1)/2 


/ U\  \lvl~ '-)/ A  _ 

P{fJ>\s,  M )  =  exp  (-(/i  +  s))  IM- 1  (2 y/jls) 
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(49) 


where  the  Bessel  function  is  expanded  using  the  series  [Gradshteyn  &  Ryzhik  1994,  §8.445]: 


!l{z)  = 
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/c=0 
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2;  \  L+2k 


(50) 


If  the  SIR  for  a  constant  target  is  given  by 
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M(A2) 

X+Pn 


(51) 


where  (■)  denotes  the  expectation.  If  K  denotes  the  fluctuation  parameter  of  a  chi-squared 
distributed  target  presented  below,  then  the  probability  of  detection  for  a  constant  target,  P^k 
(where  K  — >  oo)  can  be  given  in  terms  of  the  normalised  threshold  t  with  s  =  S, 
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(52) 


where  the  relationship  in  Equation  48  was  used  to  evaluate  the  integral  and  the  manipulation  of 
the  series  makes  use  of  the  identity 
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(53) 
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The  series  in  Equation  52  is  difficult  to  evaluate,  but  effective  methods  were  devised  by  Shnidman 
[Shnidman  1976,  Shnidman  1989].  Pd.oo  can  also  be  expressed  in  terms  of  the  generalised  Marcum 
Q-function  [Shnidman  1976], 


where 


QM(m,n ) 


The  solution  for  the  Pfa  in  Equation  46  can  be  recovered  by  setting  S  =  0  in  Equation  52.  For 
most  radar  targets  of  interest  however,  the  target  reflection  will  not  be  constant.  The  following 
section  looks  at  a  generalisation  of  this  model. 


3.2  Detection  with  a  fluctuating  target 


A  target  will  typically  fluctuate  over  the  integration  period  and  from  scan  to  scan  (different 
integration  periods).  Shnidman  [1995]  describes  four  different  fluctuation  models  including 
chi-square,  log-normal,  Weibull  and  Ricean  distributions.  The  chi-square  distribution  is  the  most 
popular  however  as  it  can  be  used  to  represent  the  well  known  Swelling  models.  With  this  model, 
the  target  fluctuation  with  parameter  K  is  given  by 

F(s|S'K>  =  f(7o(§)  exp(“^)'  <56) 

To  include  this  in  the  probability  of  detection  calculation,  it  is  multiplied  with  the  target  plus  noise 
PDF  and  then  the  product  is  integrated  over  the  random  variable  s, 

/*oo  /*oo 

PdtK(r\S,M)  =  J  P(s\S,  K)  J  P(p\s,  M)d)j,ds 

P{s\S,K)Pd,00(T\s,M)ds.  (57) 

A  general  solution  is  given  by  Shnidman  [1995], 
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(58) 


where  the  Pfa  is  identical  to  the  non-fluctuating  target  case.  This  series  is  also  difficult  to  evaluate, 
but  an  effective  method  is  given  by  Shnidman  [1995]. 

There  are  also  a  number  of  special  cases  depending  on  the  choice  of  the  parameter  K.  For  example, 
when  K  — y  oo  this  equation  reverts  back  to  the  constant  target  case  from  the  previous  section. 
Other  studies  have  looked  at  the  region  where  0  <  K  <  1,  where  the  models  are  known  as 
Weinstock  fluctuations,  [Weinstock  1964].  These  represent  slowly  fluctuating  targets  and  have 
been  shown  to  accurately  model  the  fluctuation  from  objects  which  are  cylindrical  in  shape.  Four 
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other  special  cases  of  this  model  are  when  K  =  1,  M,  2  and  2 M.  These  are  known  as  the 
Swerling  models  1-4  and  were  originally  proposed  by  [Swerling  1960]  and  detailed  in  the  books 
by  [Meyer  &  Mayer  1973,  Briggs  2004],  These  models  are  able  to  capture  the  majority  of  target 
RCS  fluctuations  and  are  popular  due  to  their  physical  interpretation,  either  representing  the  target 
fluctuation  pulse  to  pulse  or  scan  to  scan,  [Allen  &  Urkowitz  1993]. 

Examples  of  the  four  Swerling  models  are  shown  in  Figure  8  for  three  values  of  the  SIR.  Cases  1 
and  3  are  clearly  more  spread  out  due  to  the  larger  number  of  scatterers  being  observed  over  the 
slower  observation  time.  Conversely,  cases  2  and  4  have  narrower  distributions  due  to  the  faster 
change  of  RCS.  As  the  SIR  increases,  the  distribution  becomes  more  spread  and  shifts  to  the  right. 


Figure  8:  Swerling  model  comparison  for  three  different  values  of  SIR:  ( — )  case  case  2, 

(-  -)  case  3,  ( — )  case  4,  top:  S  =  —  5  dB,  middle:  S  =  0  dB,  bottom:  5  =  5  dB. 


Swerling  model  -  case  1,  (A  =  1) 


The  first  two  Swerling  cases  relate  to  targets  with  an  exponential  fluctuation.  They  represent 
several  independently  fluctuating  reflectors  of  approximately  equal  area.  This  would  apply  to 
objects  which  are  large  compared  to  the  radar  wavelength  (and  not  shaped  too  much  like  a  sphere), 
[Swerling  I960].  For  the  first  case,  a  constant  RCS  is  assumed  during  a  single  scan  (group  of 
pulses)  with  changes  from  scan  to  scan.  It  is  suitable  for  modelling  targets  viewed  with  a  high 
PRF  and  scan  rate.  The  If  solution  is  given  by, 


Pdp(r\S,M)  = 


T(M  —  1,  r) 


r(M  - 1) 

x  |  i  -  r  i  m  -  i 


+  1  +  7; 
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'  1  +  (1/S) 
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(59) 
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Swerling  model  -  case  2,  (K  =  M) 

The  Swerling  case  2  also  assumes  the  targets  have  an  exponential  fluctuation,  but  now  with 
independent  fluctuations  from  pulse  to  pulse.  It  is  suited  to  scenarios  where  targets  are  viewed 
with  a  low  PRF  or  when  small  changes  in  orientation  cause  significant  changes  in  the  illuminated 
area.  The  latter  can  arise  from  reflections  from  a  long  thin  object  subjected  to  a  high  frequency 
signal,  [Meyer  &  Mayer  1973].  For  this  case,  the  Pj  solution  is  given  by, 

PdM(r\S,  M)  =  1  -  r  (m,  1+ts/m )  /r (M).  (60) 


Swerling  model  -  case  3,  (A  =  2) 


The  target  models  defined  by  cases  3  and  4  represent  a  chi-square  distribution  with  4  degrees  of 
freedom  and  hence  there  is  not  as  much  randomness  as  the  first  two  cases.  The  model  effectively 
represents  one  large  reflector  combined  with  a  number  of  small  reflectors  or  just  one  large  reflector 
subject  to  small  changes  in  orientation.  Case  3  represents  a  target  with  constant  RCS  during  a 
single  scan  with  changes  from  scan  to  scan.  The  I\\  solution  is  given  by, 
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(61) 


Swerling  model  -  case  4,  (K  =  2 M) 


This  case  is  similar  to  the  previous  but  now  represents  a  target  with  independent  fluctuations  from 
pulse  to  pulse.  The  Pj  solution  is  given  by, 


where 
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3.3  Detection  in  correlated  clutter 

In  many  radar  systems,  there  will  be  correlation  between  pulses.  There  have  been  a  number 
of  different  attempts  to  model  the  correlation  of  both  the  clutter  and  the  target,  but  by  far  the 
most  useful  is  the  method  by  Ward,  Tough  &  Watts  [2013,  Section  12.4]  (WTW).  The  technique 
accounts  for  the  correlation  by  defining  1  <  L  <  M  as  a  measure  of  the  effective  number  of 
independent  clutter  samples.  L  is  not  necessarily  an  integer,  and  ranges  from  1  for  completely 
correlated  clutter  up  to  the  number  of  pulses  M  for  completely  uncorrelated  clutter.  Section  3.3.1 
describes  two  definitions  of  L  which  have  been  used  in  the  literature.  Section  3.3.2  then  presents 
the  full  derivation  of  the  WTW  detection  method  with  implementation  details  in  Section  3.3.3. 
The  final  Section  3.3.4  then  shows  how  the  WTW  detection  method  can  be  supplemented  with  a 
different  calculation  method  when  the  CNR  is  high. 


3.3.1  Effective  number  of  looks 


The  WTW  detection  method  requires  a  way  to  relate  the  ACF  to  an  ‘effective  number  of  looks’. 
There  are  two  alternative  methods  which  have  been  proposed  in  the  literature.  The  first  is  by 
Barton  [1969]  who  just  considered  the  ratio  of  the  correlation  time  to  the  observation  time.  Kanter 
[1986]  derived  the  effective  number  of  looks  from  this  ratio  in  terms  of  the  correlation  value,  p, 

L'  =  min  [M,  1  -  (AT  -  1)  ln(p)] .  (64) 


where  0  <  p  <  1  is  a  single  parameter  used  to  relate  the  de-correlation  time  model  to  a  correlation 
value.  This  can  be  derived  from  a  Gaussian  ACF.  Consider  the  discretised  ACF  for  the  nth  pulse, 
where  the  sample  rate  is  the  PRF,  /prf 


Pn  =  exp 


(65) 


The  second  method  is  the  ratio  of  the  squared  mean  to  the  variance  for  correlated  gamma  random 
variables  and  relates  to  the  ACF  directly.  The  result  was  originally  obtained  by  Kotz  &  Neumann 
[1963]  and  later  derived  independently  by  Kanter  [1986].  If  the  values  of  the  ACF  are  represented 
by  pn,  the  effective  number  of  looks  L  <  M  is  given  by 


M2 

M  +  2En=i\M-n)\pn\2 


(66) 


This  expression  has  a  better  statistical  justification  and  is  more  general  since  all  the  correlation  lags 
are  used.  For  comparison.  Figure  9  shows  the  number  of  effective  looks  for  each  case  using  the 
Gaussian  ACF.  Although  both  offer  similar  results,  the  second  expression  has  a  smoother  response 
as  the  correlation  increases. 
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L'  looks  L  looks 


0  0.2  0.4  0.6  0.8  1  0  0.2  0.4  0.6  0.8  1 

Correlation  Correlation 


Figure  9:  Effective  number  of  looks  comparison  with  varying  correlation  level  and  number  of 
looks,  M:  left  -  iJ ,  right  -  L. 


3.3.2  Ward,  Tough  and  Watts  calculation  method 

The  WTW  calculation  method  is  described  in  [Ward,  Tough  &  Watts  2013,  Section  12.4]  and 
provides  a  way  of  calculating  the  detection  probability  for  both  uncorrelated  and  correlated  clutter 
with  independent  noise  and  target  returns.  The  algorithm  introduces  the  quantity  (3  as  the  target 
plus  clutter  summed  over  L  independent  looks.  If  the  SIR  parameters  s  and  S  are  now  normalised 
by  only  the  noise  power  pn 


s  = 
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mi 


s  = 


M  {A2) 


Pn 


and  the  speckle  mean  power,  x  is  re-normalised  according  to 


(67) 

(68) 
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(69) 


then  the  probability  density 
2013,  Equation  12.33] 

P(P\S,<*,L) 


function  of  j3  is  an  L  pulse  Rice  distribution  [Ward,  Tough  &  Watts 
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(70) 


where  the  Bessel  function  was  expanded  using  Equation  50.  With  this  normalisation,  the  target, 
clutter  speckle  and  noise  returns  are  given  by 


M 


m=  1 


(71) 
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with  the  PDF  given  by  [Ward,  Tough  &  Watts  2013,  Equation  12.36]: 
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(72) 


Once  the  fluctuations  due  to  /3  and  s  are  included,  the  PDF  then  becomes 

/•oo  poo 

P(/i|M,L,a,5,iT)  =  /  /  P(A|/3,M)P(/3|s,a,L)P(s|5,iT)d)5ds  (73) 

7o  7o 

where  the  target  fluctuation,  P(s\S,K)  is  given  in  Equation  56.  The  solution  of  this  PDF  for 
the  cases  of  correlated  clutter  and  noise  and  a  target  present  in  correlated  clutter  can  be  found  in 
Appendix  A.  It  is  also  not  immediately  obvious  that  this  PDF  is  equivalent  to  the  uncorrelated 
case  when  the  effective  number  of  looks  L  equals  the  number  of  pulses  M.  Appendix  B  presents 
a  consistency  check  for  a  non-fluctuating  target  which  verifies  that  these  two  PDFs  are  in  fact 
equivalent. 

Calculating  of  the  probability  of  detection  therefore  requires  the  following  integral  to  be  evaluated 
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(74) 


where  the  threshold,  f  is  now  normalised  by  only  the  noise  power  and  does  not  depend  on  the 
mean  clutter  power  x. 

f  =  Y/pn.  (75) 


For  the  general  fluctuating  target  case,  the  expectation  in  Equation  74  can  be  expanded  as  [Ward, 
Tough  &  Watts  2013,  Equation  12.38  and  12.39] 
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where  the  definite  integral  [Gradshteyn  &  Ryzhik  1994,  §3.351/3]  has  been  used: 
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The  second  integral  can  then  be  evaluated  as 
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where  the  hypergeometric  function  2  Ei  is  obtained  from  [Gradshteyn  &  Ryzhik  1994,  §9.14  and 
§9.131/1]: 
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with  7]  =  — - and  (•)&  is  the  Pochhammer  symbol,  (x)fc  =  T(x  +  k)/T(x).  With 

(S  +  aK)(a  +  1) 

this  manipulation.  Equation  78  can  be  written  as 
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If  the  recurrence  relation  for  the  hypergeometric  function  [Gradshteyn  &  Ryzhik  1994,  §9.137/3] 
is  compared  with  the  recurrence  relation  for  the  function  T  [Ward,  Tough  &  Watts  2013,  Equation 
12.40],  it  can  be  established  that  \P( m ,  L,  I\,  z)  =  2-Pt(iT,  —m,  L,  —z). 


For  the  constant  target  case,  P(f3\S,a,  L)  is  obtained  from  Equation  70  with  s  =  S  and  the 
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expectation  in  Equation  74  becomes 
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(81) 


where  \F\  is  the  confluent  hypergeometric  function  and  the  relationship  [Gradshteyn  &  Ryzhik 
1994,  §9.212/1]  has  been  used: 


\F\{a,r)\  z)  =  e2iFi(7  —  a,  7;  —z). 


(82) 


If  the  recurrence  relation  for  the  function  <f>  [Ward,  Tough  &  Watts  2013,  Equation  12.42]  is 
compared  with  the  recurrence  relation  for  the  confluent  hypergeometric  function  [Gradshteyn  & 
Ryzhik  1994,  §9.212/4],  it  can  be  established  that  4>(m,  L,  z)  =  \Fi(—m,  L ;  —z). 

This  result  can  also  be  expressed  in  terms  of  a  Laguerre  polynomial  by  using  the  relationship 
[Gradshteyn  &  Ryzhik  1994,  §8.972/1]: 


m  +  n 


m 


iFi(—m,  7i+l;  x). 


cw  = 

With  this  relationship.  Equation  81  can  be  written  as 

_e-S,{a+l)cL-l  (_ 


(83) 


am 

-f)P_ 


a 


S 


m\  /  (a  +  l)L+m 

Then  the  identity  [Gradshteyn  &  Ryzhik  1994,  §8.975/1] 


a(a  +  1) 


(84) 


X>M*m  =  (l  — AT^exp 


m= 0 


XZ 

z-  1 


can  be  used  to  verify  that 


E 

771=0 


e~P  /3r‘ 

m\ 


=  1. 


(85) 
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3.3.3  Evaluating  the  detection  probability 


Before  evaluating  the  detection  probability,  Chernoff  bounds  [Shnidman  1976,  Wozencraft  & 
Jacobs  1965]  are  first  calculated  to  check  whether  the  result  can  be  approximated  as  either  0  or  1. 
Otherwise,  the  detection  probability  in  Equation  74  can  be  evaluated  by  re-writing  it  as  the  sum 
of  the  product  of  two  series,  Am  and  Cm: 


m  i 


Pd(?\(3,M)=  Y  AmC„ 


(87) 


m=mo 


where  mo  and  mi  are  the  lower  and  upper  summation  limits.  Details  of  both  the  Chernoff  bounds 
and  the  upper  and  lower  limits  are  given  in  Appendix  C.  The  first  series  in  the  sum  is  given  by 


Amo  =  exp  (mo  Inf  -  f  -  lnT(mo  +  1)) , 
Am+ 1  —  f  Amj  (m  T  1) 


(88) 


and  does  not  need  to  be  recalculated  for  different  values  of  x  since  the  threshold  f  is  independent 
of  the  clutter  power  x.  The  second  series  is 


1 An .  —  < 


L 


m  <  M 


m—M 


1  —  Bk.  m  >  M 


(89) 


k=0 


where  Bk  denotes  the  expectation 


Bu  = 


k\ 


(90) 


The  coefficients  for  the  two  target  cases  in  Equations  80  and  8 1  can  be  evaluated  using  recurrence 
relations  [Ward,  Tough  &  Watts  2013],  For  the  fluctuating  target  in  Equation  80, 


B0  = 


1 
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(91) 
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.  For  the  constant  target  case  in  Equation  8 1 , 
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where  Z  =  — - Even  with  high  precision  arithmetic,  the  series  in  Equations  91  and  92 

a(a  +  1) 

are  prone  to  underflow  for  large  S/a.  To  avoid  this,  a  method  suggested  in  Section  IV-A  of 
[Shnidman  1995]  is  used.  First  ln(i?o)  is  calculated  separately  and  stored  in  a  variable  no-  The 
series  is  started  at  1  and  each  term  is  checked  to  see  if  it  exceeds  e25.  If  it  does,  all  terms  are 
divided  by  e25  and  25  is  added  to  vo-  Once  all  required  terms  have  been  computed,  the  sum  of  the 
series  is  multiplied  by  exp(no).  In  the  constant  target  case  the  recurrence  for  Bk  is  the  <7 -function, 
and  in  the  fluctuating  target  case  the  /i-function,  both  defined  by  Shnidman  [1995].  However  the 
application  of  these  functions  is  different  here.  The  clutter  model  of  [Shnidman  1995]  is  compared 
with  the  compound  K  distribution  in  [Shnidman  2005], 

Note  that  this  calculation  method  can  also  be  used  to  calculate  the  Tfa  by  setting  S  =  0  and  using 
the  coefficients: 


B0 

Bk+i 


1 

(a  +  l)L  ’ 

L  +  k  a 

- -DL. 

k  +  1  a  +  1 


(93) 


3.3.4  Detection  at  high  CNR 


As  the  CNR  increases,  the  method  of  Ward,  Tough  &  Watts  [2013]  becomes  computationally 
intensive  and  fails  completely  when  there  is  only  clutter  present.  The  method  of  [Shnidman  1976, 
Shnidman  1989,  Shnidman  1991]  can  be  used  if  it  is  modified  for  a  real  effective  number  of 
looks  L  >  1  in  place  of  an  integer  number  M.  Consider  a  constant  target  where  the  target  plus 
interference  PDF  is  described  by  Equation  49.  The  probability  of  detection  is  then  given  by 
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(94) 


If  L  =  l  +  6,  where  l  =  floor (L),  then  using  an  identity  for  the  incomplete  gamma  function 
[Gradshteyn  &  Ryzhik  1994,  §8.356/5]: 
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Then  by  making  the  substitution  m  =  5  +  n  +  1,  where  the  index  m  is  a  real  number,  gives 


Pd(r\x,L,5)  = 


e~ssk  (  T(6  +  1,t) 


L  +  k—1  _rr 


k\ 


+  E 


e  '  r" 


k= 0 

00  -,—s  ok 


r^  +  1)  JsUnm  +  l), 


(T(6  +  1,t) 

^  k\ 


L—l 


e-rTm  L+k ^  1  e~TTm 


k= 0 

F(5  +  l,r) 


T(5  +  l)  +  2  T(m+1)+  2  r(m  +  l) 

v  '  m=S+ 1  v  '  m=L  y  , 


L—l 


Ee  / 

T(m  +  1)  '  ^  fc! 

x  '  m=5+l  K  ’  k= 0 


00  p—s„k  L+k  1  —r  m 

EC  o  C  / 

-TT-  E 


m=L 


L—l 

+  E 


T{S  +  1,t) 

r(<i  +  i)  1  r(m  +  i)  1  ^r(m  +  i) 

v  y  m=<5+l  v  7  fc=0  v  y 


+  E 


T(m  +  1) 

'  m—L  _g  fc' 

!-  E e  8 


A:=0 


fc! 


(96) 


This  expression  differs  from  the  original  summation  in  Equation  52  in  two  respects:  the  first  term 
is  new  and  the  index  m  is  a  real  number  rather  than  an  integer,  beginning  at  m  =  <5  +  1  and 
incremented  by  1  for  each  successive  term  in  the  sum.  When  L  is  an  integer,  5  =  0  and  the  first 
term  is  just  e~r  as  in  the  original  series.  For  large  r  the  original  sum  must  be  started  at  some  term 
no  >  0  to  avoid  underflow  and  loss  of  significance  (see  Appendix  C.2  for  details).  In  this  case  the 
first  term  in  Equation  96  is  negligible  and  the  sum  is  started  at  mo  =  no  +  5. 

The  same  modification  can  easily  be  applied  to  the  series  in  Equation  58  for  a  fluctuating  target. 
Note  that  the  target  fluctuation  parameter  K  is  the  same  as  for  detection  in  uncorrelated  clutter,  so 
for  example  K  =  M  for  a  Swelling  2  target  in  both  cases.  Figure  10  shows  the  effect  of  varying 
the  clutter  to  noise  ratio  on  the  probability  of  detection  for  a  Swerling  2  target  in  exponentially 
distributed  clutter  with  M  =  10  and  L  =  1.5.  The  dashed  black  line  is  calculated  using  the 
method  described  here  for  clutter  only  (infinite  CNR)  and  the  other  curves  are  calculated  with  the 
method  of  Ward,  Tough  and  Watts. 

To  determine  the  optimal  CNR  for  switching  between  methods,  the  maximum  SIR  error  between 
the  WTW  and  the  modified  Shnidman  method  has  been  measured.  As  described  in  the  following 
section,  both  of  these  calculation  methods  can  be  extended  to  compound  distributions.  The  Pareto 
distribution  has  been  used  to  model  the  clutter  with  a  shape  of  3  and  the  number  of  looks,  M  =  5. 
To  establish  the  worst  error,  the  speckle  is  totally  correlated  with  p  =  1  and  hence  L  =  1.  Fig. 
1 1  shows  the  maximum  SIR  error  and  the  mean  run  time  averaged  over  3  runs  as  a  function  of 
the  CNR.  As  expected,  the  SIR  error  is  0.7  dB  when  the  CNR  is  at  -10  dB  due  to  the  incorrect 
modelling  of  the  correlated  thermal  noise.  However,  this  quickly  reduces  and  becomes  less  than 
0. 1  dB  at  0  dB.  While  the  WTW  detection  method  is  more  accurate,  it  is  far  more  computationally 
intensive.  At  the  same  transition  point  when  the  CNR  is  0  dB,  the  run  time  is  around  10  s  and 
increases  to  70  s  at  10  dB.  For  comparison,  the  mean  run  time  for  the  modified  Shnidman  algorithm 
with  any  value  of  CNR  is  0.14  s.  From  this  result,  a  realistic  transition  point  between  algorithms 
would  be  at  0  dB  or  when  the  clutter  starts  to  dominate  the  interference. 
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Figure  10:  Pci  with  different  clutter  to  noise  ratios  for  a  Swerling  2  target  in  exponentially 
distributed  clutter  and  noise  after  non-coherent  integration  of  10  pulses  with  a  Pfa 
of  10  K  The  effective  number  of  independent  clutter  samples  is  1.5. 


Figure  11:  Maximum  SIR  error  between  detection  algorithms  (top)  and  mean  run  time  for  WTW 
algorithm  (bottom). 
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3.4  Detection  in  compound  clutter 

Thus  far  the  clutter  power  x  has  been  treated  as  a  constant.  In  the  compound  model  x  is  a  random 
variable  with  a  PDF  Px(x).  The  global  probability  of  detection  P'A  is  obtained  by  integrating  the 
local  PA  over  this  PDF.  The  effect  of  this  ‘smoothing’  is  to  degrade  the  non-coherent  integration 
by  up  to  20  dB  when  compared  to  homogeneous  Gaussian  clutter  [Allen  &  Urkowitz  1993]. 

With  this  extension,  the  local  SIR,  S  is  now  replaced  by  a  global  SIR,  S'  and  the  local  threshold, 
r  by  a  global  threshold,  t', 


ff  M(A*)' 

Pc  +  Pn' 

,  Y 

T  =  - 

Pc+Pn 


(97) 

(98) 


These  can  both  be  rewritten  in  terms  of  the  short  scale  definitions, 

(99) 

(100) 

The  Pd  solution  is  determined  by  multiplying  the  PDF  of  the  mean  speckle  power,  Px  (x)  and  then 
integrating  over  x. 


S' (Pc+Pn) 
X+Pn 
r’jpc  +Pn) 
X+Pn 


P^k(t'\S',M)  =  /  Px(x)PdiK 


(  t' (pc  +  Pn)  \S' (pc  +  Pn) 
V  X+pn  X+Pn 


M  dx. 


(101) 


To  the  authors’  knowledge,  there  is  no  documented  analytic  solution  to  this  integral  and  numerical 
integration  must  be  used  for  evaluation  as  described  in  Section  2.3  of  [Bocquet  2012].  Different 
target  fluctuations  can  be  modelled  by  changing  the  value  of  K  as  in  the  previous  section.  The 
solution  for  PA  can  be  expressed  in  terms  of  t',S'  and  the  CNR,  C' .  If  the  integration  variable  is 
redefined  as  u  =  x/pc. 


+^S',M,C)  =  l  -P.  dx 

V(C"  +  1).S,(C"  +  1) 


Pu(u)PdtK 


t'(C'  +  1)  ,  S'(C'  +  1) 

'cC  +  1 
,  M)  du 


uC'  +  l  1  uC'  +  l 


where  Px(u )  =  1  /pcPx(x/pc)  and  dx  =  pcdu. 

The  probability  of  false  alarm  can  be  calculated  similarly, 

1 


PUIr'IM,  C)  =  jT  W  (m,  TgAT 


du 


(102) 


(103) 


where  again  numerical  integration  can  be  used  for  evaluation. 

For  the  WTW  method,  the  same  procedure  is  applied,  except  that  the  parameters  S  and  t  are 
normalised  by  the  noise  power  alone  and  hence  are  independent  of  x.  If  these  parameters  are 
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renormalised  by  the  sum  of  clutter  and  noise  powers  we  obtain 
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C  +  1’ 
f 

c>  + 1' 


The  dependence  on  x  is  only  through  the  parameter  a, 

Mx  M 
Lpn  L 


uC'. 
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4  Detection  performance  analysis 

This  section  presents  target  detection  performance  results  using  the  empirical  models  from  Section 
2.  In  these  results,  the  Pareto  distribution  acts  as  a  proxy  for  the  data  as  it  has  been  shown  to 
match  extremely  well,  [Rosenberg  &  Bocquet  201 3«] .  By  accounting  for  temporal  correlation 
and  comparing  the  differences  in  detection  performance  between  the  K  and  Pareto  distributions, 
this  analysis  also  highlights  the  importance  of  using  the  correct  model  when  assessing  detection 
performance.  The  outcome  of  this  section  is  an  accurate  description  of  the  detection  performance 
given  a  variety  of  target  fluctuation  models,  collection  geometries  and  environmental  conditions. 

As  a  basis  for  the  modelling  results,  a  common  set  of  radar  parameters  are  described  in  Table 
6.  The  modelled  outputs  described  in  Section  2.3.2  are  then  shown  in  Tables  7-9  for  a  variety 
of  geometries,  polarisations  and  sea-states.  The  trends  for  the  shape  parameters  include  a  linear 
increase  with  grazing  angle,  an  increase  with  sea-state  and  a  maximum  in  the  upwind  direction. 
Also,  due  to  the  effect  of  the  platform  motion,  the  overall  decorrelation  time,  T' ,  from  Equation  29 
is  typically  smaller  than  the  temporal  decorrelation,  T,  from  Equation  30.  The  clutter  correlation 
is  determined  using  the  overall  decorrelation  time  with  the  Gaussian  model  in  Equation  65  and  the 
effective  number  of  looks  is  determined  with  Equation  66. 

The  first  section  4. 1  shows  the  false  alarm  rate  for  the  K  and  Pareto  distributions.  Section  4.2  then 
describes  the  Pd  results  for  a  selection  of  fluctuating  target  models,  false  alarm  rates,  geometry 
and  sea-states.  This  then  leads  into  the  final  Section  4.3  which  looks  at  the  minimum  detectable 
target  RCS. 

Table  6:  Common  simulation  input  parameters 


Parameters 

Symbol 

Value 

Centre  frequency 

fc 

10  GHz 

Bandwidth 

B 

200  MHz 

Pulse  width 

TP 

5  /is 

PRF 

/PRF 

500  Hz 

Peak  transmitted  power 

Pt 

20  kW 

Antenna  gain  (one  way) 

G 

35  dB 

Combined  losses 

La  +  Lg 

7  dB 

Noise  figure 

Fn 

4  dB 

Reference  temperature 

To 

290  K 

Azimuth  beam  width  two-way  3  dB 

03dB 

1° 

Doppler  spectrum  mean 

fo 

0  Hz 

Platform  velocity 

Vp 

100  m/s 

Slant  range 

R 

20  km 

Elevation  /  grazing  angle 

9 

30° 

Azimuth  angle 

0 

0° 

Wind  swell  angle 

0° 

Sea-state 

5 

3 

Number  of  looks 

M 

10 

Probability  of  false  alarm 

Ufa 

10“5 

Probability  of  detection 

Pd 

0.5 

Polarisation 

HH 

Target  fluctuation  model 

Swerling  case  2 
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Table  7:  Calculated  model  parameters  for  varying  azimuth  angle  with  30°  grazing,  sea-state  3. 


Polarisation 

Azimuth  angle 

CNR  (dB) 

V 

a 

T  (ms) 

T'  (ms) 

HH 

Upwind 

14.60 

2.62 

4.10 

10.11 

4.66 

Crosswind 

9.23 

1.51 

2.86 

Downwind 

10.68 

3.29 

4.53 

HV 

Upwind 

3.22 

7.18 

8.15 

11.43 

4.78 

Crosswind 

0.63 

1.36 

2.53 

Downwind 

2.68 

5.62 

7.91 

Upwind 

20.92 

31.77 

34.94 

VV 

Crosswind 

15.73 

16.52 

16.55 

12.75 

4.86 

Downwind 

19.79 

36.07 

37.12 

Table  8:  Calculated  model  parameters  for  varying  grazing  angle  with  sea-state  3,  upwind 
direction. 


Polarisation 

Grazing  angle 

CNR  (dB) 

V 

a 

T  (ms) 

T'  (ms) 

HH 

to 

o 

o 

10.16 

0.60 

2.37 

10.11 

4.37 

o 

o 

co 

14.60 

2.62 

4.10 

4.66 

45° 

21.61 

11.40 

14.05 

5.43 

HV 

to 

o 

o 

1.27 

3.24 

4.20 

11.43 

4.46 

o 

o 

co 

3.32 

7.18 

8.15 

4.78 

45° 

6.74 

15.88 

19.20 

5.61 

to 

o 

o 

18.40 

30.60 

33.70 

4.53 

VV 

CO 

o 

o 

20.92 

31.77 

34.94 

12.75 

4.86 

45° 

25.05 

32.99 

36.23 

5.75 

Table  9:  Calculated  model  parameters  for  varying  sea  state  with  30°  grazing,  upwind  direction. 


Polarisation 

Sea  state 

CNR  (dB) 

V 

a 

T  (ms) 

T'  (ms) 

HH 

1 

5.96 

0.42 

2.37 

14.22 

4.98 

3 

14.6 

2.62 

4.10 

10.11 

4.66 

6 

20.1 

4.13 

5.15 

3.65 

3.00 

HV 

1 

-6.40 

0.37 

3.05 

14.94 

4.96 

3 

3.32 

7.18 

8.15 

11.43 

4.78 

6 

9.45 

10.82 

10.25 

5.74 

3.88 

VV 

1 

11.32 

22.07 

23.06 

15.63 

4.98 

3 

20.92 

31.77 

34.94 

12.75 

4.86 

6 

26.98 

29.92 

32.17 

7.83 

4.36 

4.1  Probability  of  false  alarm 

The  probability  of  false  alarm  is  related  to  the  CCDF  as  shown  in  Equation  39.  In  a  detection 
scheme,  the  desired  Pfa  is  specified  and  the  threshold  must  be  found  by  inverting  this  equation. 
Using  the  parameters  above,  Figure  12  shows  a  comparison  of  Pfa  for  the  K  and  Pareto 
distributions.  A  range  of  looks,  M  (non-coherent  sum  of  pulses)  has  also  been  shown  to 
demonstrate  the  effect  of  correlation. 
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There  is  a  clearly  a  big  difference  in  the  HH  and  HV  results  with  the  K-distribution  having  a 
lower  threshold  by  approximately  5  dB  and  2  dB  respectively  at  10“5,  while  there  is  only  a  minor 
difference  in  the  VV  channel.  The  effect  of  non-coherently  integrating  pulses  is  to  reduce  the 
threshold.  This  effect  is  more  pronounced  in  the  VV  polarisation  where  the  threshold  is  reduced 
by  6  dB  with  50  pulses  integrated. 


M  =  50 


Normalised  threshold  (dB) 


Figure  12:  Pfa  comparison  with  number  of  looks:  ( — )  K,  (-  -)  Pareto. 


4.2  Probability  of  detection 

Once  a  Pfa  is  chosen  and  the  threshold  is  found,  the  detection  probability,  Pd  can  be  calculated  as 
a  function  of  the  Signal  to  Interference  Ratio  (SIR).  The  next  set  of  results  in  Figures  13-17  show 
variations  in  the  target  model,  number  of  looks,  geometry  and  sea-state. 

The  target  fluctuation  models  can  all  be  derived  from  the  chi-squared  distribution  model  in 
Equation  56.  A  non-fluctuating  Marcum  target  model  is  found  when  I\  — >  oo.  Using  the  Pareto 
distribution,  the  Pd  result  is  shown  along  the  top  row  of  Figure  13  for  M  =  1, 10  and  50  looks 
and  a  Pfa  of  10  5.  As  the  number  of  looks  increases,  the  curve  shifts  to  the  left  and  the  required 
SIR  for  a  contant  Pd  is  reduced.  The  Swerling  fluctuation  models  1  and  2  refer  to  the  cases  where 
K  =  1  and  K  =  M  and  the  target  model  becomes  an  exponential  distribution.  Similarly,  the 
Swerling  models  3  and  4  refer  to  a  chi-squared  distribution  with  K  =  2  and  K  =  2 M.  To 
account  for  target  correlation,  an  effective  value  B  <  M  can  be  determined  using  the  relationship 
in  Equation  66  with  an  appropriate  target  ACF.  For  targets  with  an  exponential  distribution,  the 
fluctuation  parameter  is  related  by  K  =  B,  while  for  chi-squared  distributions,  the  relationship  is 
K  =  2 B.  If  a  Gaussian  ACF  is  assumed  for  the  target  as  in  Equation  65,  then  the  detection 
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performance  can  be  calculated  as  a  function  of  the  target  correlation,  pt(n )  =  p"  for 
n  =  0, . . . ,  M  —  1.  For  pt  =  0,  the  target  is  uncorrelated  and  B  =  1,  while  pt  =  1  implies  total 
correlation  and  li  =  M.  Figure  13  shows  a  comparison  for  both  exponential  and  chi-squared 
target  models  as  the  number  of  looks  increases.  The  increasing  number  of  looks  has  again  shifted 
the  curves  to  the  left,  while  an  increasing  target  correlation  reduces  the  slope  of  the  curve.  In  all 
cases,  the  chi-squared  result  has  a  slightly  steeper  slope  than  the  exponential  fluctuation.  Since 
the  target  correlation  cannot  be  known  in  advance,  only  the  uncorrelated  exponential  fluctuation 
(Swerling  2)  model  is  used  for  the  remainder  of  these  comparisons.  Table  10  shows  the  minimum 
required  SIR  for  a  Pd  =  0.5. 

Using  a  model  also  allows  us  to  extrapolate  down  to  lower  Pfa  levels  than  could  be  used  in  a 
real  detection  scheme.  The  result  in  Figure  14  shows  the  K  and  Pareto  distribution  results  for  a 
Pfa  of  10-4,  1 0  r'  and  10fi.  For  both  the  FTFT  and  F1V  polarisations,  there  is  a  clear  difference 
between  the  two  distribution  models  which  increases  as  the  Pfa  reduces.  For  the  HF1  polarisation, 
the  difference  in  the  SIR  at  a  Pd  of  0.5  is  up  to  5.4  dB,  while  this  difference  is  reduced  for  the 
cross-pol  polarisation  and  almost  zero  for  VV.  Table  1 1  shows  the  required  SIR  for  each  result 
measured  at  a  Pd  of  0.5. 

The  next  three  results  show  a  comparison  with  geometry  and  sea-state.  Figure  15  first  shows  the 
variation  in  azimuth  direction.  The  biggest  variation  between  the  K  and  Pareto  distributions  is 
found  for  HF1  and  HV  cross  wind  directions  where  the  SIR  difference  is  up  to  7.4  dB.  There  is  little 
difference  in  the  VV  polarisation,  where  the  Pareto  matches  the  K  distribution  result.  Table  12 
shows  the  required  SIR  for  each  result  measured  at  a  Pd  of  0.5. 

Figure  16  shows  a  comparison  of  20°,  30°  and  45°  grazing.  The  biggest  difference  between  the 
distributions  is  found  in  the  F1H  polarisation  where  the  SIR  difference  is  up  to  5.1  dB  for  20° 
grazing.  The  mismatch  is  reduced  for  the  VV  polarisation  and  also  when  the  grazing  angle  gets 
higher.  Table  13  shows  the  required  SIR  for  each  result  measured  at  a  Pd  of  0.5. 

The  final  result  in  Figure  17  shows  the  variation  in  sea-state.  The  largest  mismatch  of  5.1  dB 
is  found  for  the  lowest  sea-state  in  the  HF1  polarisation.  The  mismatch  reduces  as  the  sea-state 
increases  and  is  very  small  for  the  VV  polarisation.  For  the  lowest  sea-state,  the  CNR  for  the 
F1V  channel  is  very  low  and  both  distributions  are  approximately  exponential  giving  the  same 
performance  results.  Table  14  shows  the  required  SIR  for  each  result  measured  at  a  Pd  of  0.5. 


Table  10:  Required  SIR  (dB)  at  a  Pj  of  0.5.  Variation  in  target  model  and  number  of  looks 
corresponding  to  Figure  13. 


M 

1 

10 

50 

Marcum 

16.66 

14.54 

13.50 

Exponential,  pt  =  1  (Swelling  1) 

18.23 

16.11 

15.07 

Exponential,  pt  =  0  (Swelling  2) 

18.23 

14.67 

13.51 

Chi-squared,  pt  =  1  (Swerling  3) 

17.40 

15.28 

14.24 

Chi-squared,  pt  =  0  (Swerling  4) 

17.40 

14.60 

13.50 
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M  =  1  M  =  10  M  =  50 


0  10  20  30  0  10  20  30  0  10  20  30 

SIR  (dB)  SIR  (dB)  SIR  (dB) 


Figure  13:  P,j  variation  with  target  model  and  number  of  looks.  Variation  with  target  correlation: 


Pt:  (~)  0,  (-)  0.5,  (-)  0.9,  (-)  1. 


P,  =10'4  P,  =  10"5  P,  =10'6 

fa  fa  fa 


SIR  (dB)  SIR  (dB)  SIR  (dB) 


Figure  14:  Pd  variation  with  polarisation  and  false  alarm  rate:  ( — )  K,  ( — )  Pareto. 
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Table  11:  Required  SIR  (dB)  at  a  Pd  of  0.5.  Variation  in  false  alarm  rate  corresponding  to  Figure 
14. 


Pol. 

l\  = 

=  10"4 

Ufa  = 

10 b 

Ufa  = 

10"b 

K 

Pareto 

K 

Pareto 

K 

Pareto 

HH 

9.37 

11.99 

10.71 

14.67 

11.82 

17.24 

HV 

6.76 

7.61 

7.88 

9.20 

8.82 

10.66 

VV 

6.10 

6.16 

7.10 

7.22 

7.94 

8.11 

Table  12:  Required  SIR  (dB)  at  a  Pd  of  0.5.  Variation  in  azimuth  angle  corresponding  to  Figure 
15. 


Pol. 

Upwind,  f  =  0° 

Crosswind,  f  =  90° 

Downwind,  f  =  180° 

K 

Pareto 

K 

Pareto 

K 

Pareto 

HH 

10.71 

14.67 

11.80 

18.16 

10.03 

13.65 

HV 

7.88 

9.20 

10.28 

17.63 

8.11 

9.15 

VV 

7.10 

7.22 

7.66 

8.28 

7.01 

7.16 

Upwind  Crosswind  Downwind 


SIR  (dB)  SIR  (dB)  SIR  (dB) 


Figure  15:  Pd  variation  with  polarisation  and  azimuth  angle:  ( — )  K,  ( — )  Pareto. 
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Grazing  =  20°  Grazing  =  30°  Grazing  =  45° 


"O 


o  ^ - * - * - ■  ^ - * - ■ - ■  ^ - * - * - 1 

0  10  20  30  0  10  20  30  0  10  20  30 

SIR  (dB)  SIR  (dB)  SIR  (dB) 


Figure  16:  F,j  variation  with  polarisation  and  grazing  angle:  ( — )  K,  ( — )  Pareto. 
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Seastate  =  3 
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Figure  17:  F,j  variation  with  polarisation  and  sea-state:  ( — )  K,  ( — )  Pareto. 
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Table  13:  Required  SIR  (dB)  at  a  Pc\  of  0.5.  Variation  in  grazing  angle  corresponding  to  Figure 
16. 


Pol. 

e  = 

20° 

0  = 

CO 

o 

o 

e  = 

=  45° 

K 

Pareto 

K 

Pareto 

K 

Pareto 

HH 

14.45 

20.89 

10.71 

14.67 

8.40 

8.91 

HV 

8.53 

12.27 

7.88 

9.20 

7.75 

7.96 

VV 

6.94 

7.07 

7.10 

7.22 

7.51 

7.61 

Table  14:  Required  SIR  (dB)  at  a  Pj  of  0.5.  Variation  in  sea-state  corresponding  to  Figure  17. 


Pol. 

5 

=  1 

S 

=  3 

S 

=  6 

K 

Pareto 

K 

Pareto 

K 

Pareto 

HH 

15.24 

20.31 

10.71 

14.67 

8.94 

12.50 

HV 

14.93 

14.93 

7.88 

9.20 

7.48 

8.69 

VV 

7.36 

7.36 

7.10 

7.22 

6.89 

7.06 

4.3  Minimum  detectable  target  RCS 

The  final  section  looks  at  the  the  minimum  detectable  point  target  RCS  for  a  given  set  of  radar 
parameters  and  a  desired  Pfa  and  If.  These  can  be  related  using  the  radar  range  equation  as  it 
applies  to  both  the  target  and  the  sea-clutter.  For  the  sea-clutter,  the  mean  power  was  given  in 
Equation  1 .  For  the  mean  target  power, 


PtG2X2cotTpB 

(4ir)3R4LaLs 


(107) 


where  at  is  the  effective  target  RCS  inclusive  of  any  propagation  effects.  To  determine  the 
minimum  detectable  target  RCS,  crt,m in,  the  relationship  for  the  SIR  can  be  used: 


s  = 


Ps 

Pc +Pn 


(108) 


If  this  relationship  is  then  written  in  terms  of  the  minimum  detectable  target  RCS,  (Jt.mm, 


n't, min 


■SO, min  I  nc  4~ 


kT0Fn(4TrYR4LaL 

PtG2X2Tp 


(109) 


where  so,min  is  the  minimum  required  SIR.  Using  the  parameters  given  above,  Figure  18  shows  the 
minimum  SIR  and  the  corresponding  minimum  target  RCS  for  each  polarisation  as  it  varies  with 
grazing  angle.  Also  shown  in  both  results  is  the  noise  limited  case  (pc  =  0)  which  is  constant  at 
3  dB  for  the  minimum  SIR  and  -14  dB  for  the  minimum  target  RCS.  For  the  lower  grazing  angles 
in  the  HH  and  HV  channels,  there  is  a  large  difference  between  the  K  and  Pareto  distributions, 
which  reduces  as  the  grazing  angle  increases.  The  two  distributions  give  almost  identical  results 
for  the  VV  channel.  In  terms  of  detectability,  the  Pareto  result  shows  that  the  HV  polarisation  has 
the  smallest  target  RCS,  then  VV  and  HH. 
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15  20  25  30  35  40  45 

Grazing  (deg) 


15  20  25  30  35  40  45 

Grazing  (deg) 


Figure  18:  Minimum  required  SIR  (left)  and  minimum  detectable  target  RCS  (right)  with  M  =  10 
looks  and  Pfa  =  10-J.  Other  parameters  given  in  Table  6:  ( — )  K,  ( — )  Pareto,  ( — ) 
noise  limited  case. 
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5  Conclusions  and  future  work 

Parametric  modelling  is  a  useful  tool  to  predict  a  radar’s  target  detection  performance  given  a 
statistical  description  of  the  environment.  This  report  summarised  the  target  detection 
performance  of  an  X-band  surface  surveillance  radar  operating  at  medium  grazing  angles  in  the 
maritime  environment.  A  number  of  different  detection  scenarios  were  explored  with  variations 
in  the  collection  geometry,  sea-state  and  polarisation.  Due  to  its  good  fit  to  the  Ingara  sea  clutter 
data,  the  Pareto  distribution  was  used  as  a  proxy  for  the  data  with  the  K-distribution  results 
showing  the  mismatch  that  would  be  expected  in  a  detection  scenario.  The  overwhelming 
conclusion  from  this  report  is  that  the  choice  of  model  is  critical  to  accurately  determine  the 
detection  performance. 

These  scenarios  required  a  number  of  parameter  models  to  relate  the  radar  and  environmental 
characteristics  to  the  appropriate  mean  backscatter  and  the  distribution  shape  values.  Appropriate 
models  were  therefore  presented  along  with  the  equations  used  to  determine  the  performance  of 
different  fluctuating  targets  in  both  uncorrelated  and  correlated  clutter.  A  full  derivation  of  the 
Ward,  Tough  and  Watts  detection  algorithm  for  correlated  clutter  was  presented  with  a 
modification  to  make  the  technique  robust  for  high  CNRs. 

The  sea-clutter  and  target  detection  models  were  then  used  to  analyse  the  probability  of  false  alarm, 
probability  of  detection  and  the  minimum  detectable  RCS  for  a  given  set  of  radar  parameters. 
The  Pfa  results  showed  a  a  big  difference  in  the  horizontal  polarisations  with  the  K  distribution 
underestimating  the  threshold  by  up  to  5  dB,  while  there  was  only  a  minor  difference  in  the  VV 
channel.  The  effect  of  non-coherently  integrating  pulses  was  to  reduce  the  threshold. 

A  comparison  of  the  target  fluctuation  models  revealed  that  the  chi-squared  distribution  model 
has  a  slightly  lower  I\\  than  the  exponential  distribution  and  that  increasing  the  target  correlation 
reduced  the  slope.  Using  a  distribution  model  also  allows  us  to  extrapolate  down  to  lower  Pfa 
levels  than  could  be  used  in  a  real  detection  scheme.  By  comparing  Pfa  levels,  it  was  found  that 
for  the  HH  polarisation,  the  difference  in  the  SIR  between  the  K  and  Pareto  distributions  was  up 
to  5.4  dB.  This  mismatch  also  got  larger  as  the  Pfa  reduced.  When  looking  at  changes  in  the 
azimuth  look  direction,  the  biggest  variation  between  the  K  and  Pareto  distributions  was  found 
for  the  horizontal  crosswind  directions  where  the  SIR  difference  was  up  to  7.7  dB.  There  was 
little  difference  in  the  VV  polarisation,  where  the  Pareto  matched  the  K-distribution  result.  For 
changes  in  grazing,  a  difference  was  observed  between  the  distribution  models  of  up  to  6.4  dB. 
The  mismatch  is  reduced  for  the  VV  polarisation  and  also  when  the  grazing  gets  higher.  The  final 
comparison  was  for  the  sea-state,  where  the  largest  mismatch  of  5.1  dB  is  found  for  the  lowest 
sea-state  in  the  HH  polarisation.  The  mismatch  reduced  as  the  sea-state  increased  and  was  very 
small  for  the  VV  polarisation. 

The  final  result  looked  at  the  minimum  detectable  target  RCS  given  a  set  of  radar  parameters.  For 
the  lower  grazing  angles  in  the  HH  and  HV  channels,  there  was  a  large  difference  between  the  K 
and  Pareto  distributions,  which  reduced  as  the  grazing  angle  increased.  In  terms  of  detectability, 
the  Pareto  result  shows  that  the  VV  polarisation  has  the  smallest  target  RCS,  then  HV  and  HH. 

These  parameter  models  have  also  been  recently  used  for  simulating  realistic  sea-clutter  with  an 
evolving  Doppler  spectrum  along  time  and  range  [Watts  2012,  Watts,  Rosenberg  &  Ritchie  2014, 
Bocquet,  Rosenberg  &  Watts  2014].  By  injecting  false  targets  into  this  simulated  data,  detection 
performance  can  be  analysed  using  adaptive  CFAR  algorithms  and  other  ‘new’  detection  schemes. 
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Appendix  A  PDFs  for  correlated  clutter 

This  appendix  derives  the  PDFs  from  the  WTW  detection  model  in  Section  3.3.2  for  the  cases  of 
correlated  clutter  and  noise  and  a  target  in  correlated  clutter  and  noise.  The  PDFs  derived  below  are 
for  clutter  speckle,  so  numerical  integration  over  the  texture  PDF  is  required  to  obtain  the  PDF  for 
compound  clutter,  as  described  in  Section  3.4.  In  order  to  obtain  accurate  results  for  both  low  and 
high  intensity,  the  integration  interval  is  separated  into  two  parts,  with  100  point  Gauss-Legendre 
quadrature  used  on  [0, 1]  and  20  point  Gauss-Laguerre  quadrature  on  [1,  oo). 


A.l  Correlated  clutter  and  noise 


The  PDF  for  correlated  clutter  and  noise  is  evaluated  as  follows.  With  no  target  present,  s  =  0  and 
Equation  70  becomes 

aL-l  -L  -0/a 

P(P\a,L)  =  P - ^ - .  (Al) 


F  (L) 


The  PDF  can  then  be  written  as 

POO 


/»oo 
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(A2) 


and  is  the  same  PDF  obtained  by  [Swerling  1970,  Equation  54]  for  a  fluctuating  target  in  noise.  It 
is  also  the  ‘noncentral  chi-square  gamma’  PDF  used  by  Shnidman  [1999]  as  a  clutter  model.  The 
confluent  hypergeometric  function  can  be  evaluated  using  a  Taylor  series  for  small  arguments  and 
an  asymptotic  series  for  large  arguments  [Pearson  2009].  Here  the  parameters  and  argument  are 
all  positive  and  real  or  integer,  so  the  function  is  also  real  and  positive  with  no  singularities.  In  this 
case,  the  asymptotic  series  can  be  simplified  to 


1P1  (L;M-,u) 


F  (N)  L-N  u  -  L)k  (1  ~  L)k  —k 

T(L)U  6  2s  k\ 

'  J  i — n 


(A3) 


The  asymptotic  series  is  used  for  u  >  100,  and  the  Taylor  series  for  u  <  100.  If  the  Matlab 
symbolic  math  toolbox  is  available,  the  function  ‘hypergeom’  could  be  used,  but  according  to 
[Pearson  2009],  this  function  is  at  least  100  times  slower  than  direct  evaluation  of  the  series.  The 
PDF  is  obtained  by  numerical  integration  of  Equation  A2  over  the  gamma  distribution  for  the 
clutter  power  x.  The  dependence  on  x  is  only  through  the  parameter  a,  per  Equation  69. 


Figure  Al  shows  the  PDF  for  M  =  10  pulses  and  different  values  of  L,  with  a  CNR  of  0  dB 
and  K  distribution  shape  parameter  u  =  0.1.  The  clutter  correlation  mainly  affects  the  tail  of  the 
distribution,  and  hence  the  probability  of  false  alarm  for  a  given  threshold. 
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Figure  Al:  PDF  for  the  sum  of  10  returns  from  correlated  K  clutter  and  noise,  with  a  CNR  ofO  elB 
and  shape  parameter  u  =  0.1,  for  differen  t  numbers  of  independen  t  clutter  samples. 


A. 2  Target  in  correlated  clutter  and  noise 

The  PDF  for  a  target  in  correlated  clutter  and  noise  is 

roo  roc 

P(pt\M,L,a,S,K)=  /  P(ji\p,M)P(p\5,a,L)P(5\S,K)dpd. 

Jo  Jo 


(A4) 


The  integral  over  /3  can  be  evaluated  using  the  definite  integral  [Gradshteyn  &  Ryzhik  1994, 
§6.633]: 
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If  the  target  is  non-fluctuating,  s  =  S,  and  the  definite  integral  from  Equation  A7  can  be  applied 
to  obtain  the  PDF 
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If  there  is  no  target  present,  S  =  0  and  Equation  A8  reduces  to  Equation  A2.  For  a  fluctuating 
target  the  PDF  is 
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Appendix  B  WTW  detection  model  consistency  check 

This  appendix  presents  a  consistency  check  for  the  WTW  detection  model  to  show  that  in  the  case 
of  a  non-fluctuating  target,  the  M  pulse  Rice  distribution  is  recovered  when  the  effective  number 
of  looks  L  equals  the  number  of  pulses  M.  This  can  be  verified  by  integrating  the  product  of  the 
PDFs  in  Equations  70  and  72  over  (5, 


P{fi\s,a,M)=  /  P(p\p,M)P(p\s,a,M)dp 
Jo 


1^2  exp  (— jj,  —  s/a) 


j  exp  Im-i  (2 Im-  1  (~~~  )  dP 


(Bl) 


where  the  integral  in  Equation  Bl  can  be  evaluated  using  [Gradshteyn  &  Ryzhik  1994,  §6.615]: 
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11  <1 

By  making  the  following  substitutions:  /  =  jF  =  e2J F  and  g  =  jG  =  e2JG  and  noting  that 
Im(z)  =  e_m2 i  Jm  (ei-’z)  [Gradshteyn  &  Ryzhik  1994,  §8.406/1],  and  Im(—z)  =  Im  (e71"-7 z)  = 
e“ TJ Im(z)  from  the  series  [Gradshteyn  &  Ryzhik  1994,  §8.445],  Equation  B2  can  be  written  as 
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The  Rice  distribution  for  the  sum  //  of  M  radar  returns  from  a  target  in  noise  and  uncorrelated 
clutter  speckle  was  given  in  Equation  49  as 

/n\(M- 1)/2 

P(n\s,M)  =  (Z)  e~^+shM-i(2  y/JFs)  (B5) 

and  since  a  =  x/pn  when  L  =  M,  we  have  p  =  p/[a  +  1)  and  s  =  s/(a  +  1).  Then  by 
comparing  Equation  B4  with  Equation  B5,  it  can  be  seen  that  P(g\s,  a,  M)dg  =  P{p\s,  M)dp. 
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Appendix  C  Calculation  bounds  for  the  WTW  detection 

probability 

This  appendix  supplements  the  WTW  algorithm  description  in  Sections  3.3.2  and  3.3.3.  Section 
C.l  describes  the  Chernoff  bounds  which  are  used  to  detect  cases  when  the  Pd  is  less  than  the 
minimum  value  that  can  be  represented  on  a  computer  or  when  1  —  P&  is  less  than  the  smallest 
number  which  can  be  represented  with  double  precision  and  the  best  solutions  will  be  either  0  or 
1  respectively.  Section  C.2  then  derives  practical  lower  and  upper  limits  for  evaluating  the  h\\  in 
Equation  87. 


C.l  Chernoff  bounds 

Chernoff  bounds  provide  a  means  to  avoid  unnecessary  calculations  when  the  probability  of 
detection  is  close  to  zero  or  one  [Chernoff  1952,  Wozencraft  & 
Jacobs  1965,  Shnidman  1976,  Shnidman  1995,  Bocquet  2012].  They  are  a  function  of  the 
number  of  looks,  M,  the  SIR,  S,  the  threshold,  f  and  a  parameter,  A.  If  an  upper  bound  is 
desired,  then  the  value  of  A  must  be  found  which  most  closely  matches  the  inequality 

Pd(M,S,f)  <  Cb(M,S,t ,  A)  (Cl) 


or  for  a  lower  bound 


Pd(M,S,f)  >  1  —  CB(M,S,f,X).  (C2) 

To  use  the  Chernoff  bound,  the  optimal  Ao  is  calculated  for  the  appropriate  target  model  and  two 
cases  are  investigated: 

Ao  is  positive  implies  Pd  is  less  than  Cb{M,  S,  f ,  Ao)  and  we  set  Pd  to  0. 

Ao  is  negative  implies  Pd  is  greater  than  Cb(M,  S.  t,  Ao),  and  we  set  Pd  to  1. 

When  applied  to  the  compound  distributions,  the  Chernoff  bounds  are  calculated  prior  to 
numerical  integration  for  each  value  of  a  used  in  the  integration.  In  Matlab,  this  can  be  done 
simultaneously  by  treating  the  variable  a  as  a  vector  formed  from  all  the  required  values.  Even 
though  the  calculation  of  Ao  is  rather  complicated,  evaluation  of  the  Chernoff  bound  is  generally 
much  faster  than  summing  the  series  for  the  probability  of  detection,  especially  for  large  f  where 
many  terms  are  needed  for  the  series  to  converge. 


Non-fluctuating  target  in  uncorrelated  clutter 


For  a  non-fluctuating  target  in  noise,  the  Chernoff  bound  is  [Shnidman  1995,  Equation  73] 


CB(M,S,f,\ ) 


exp(— At  +  <SA/(1  —  A)) 
(1  -  \)M 


(C3) 
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Non-fluctuating  target  in  correlated  clutter 


The  bound  for  an  effective  number  of  looks  L  in  correlated  clutter  plus  noise  is  obtained  by  setting 
S  =  (3  and  integrating  over  the  probability  P(B\S,  a,  L ): 

poo 

Cb{M,  L,  a,  S,  f ,  A)  =  /  CB{M,/3,f,\)P(/3\S,a,L)dp 

Jo 

exp(— Af  +  /?A/(1  —  A))  1  / (3  A  s 


L 


(1  -  A)M 

exp  Af  —  S/a'j 


a  \SJ  \  a  J  \  a 


POO 

/  exp(-/3(A/(A  -  1)  +  1/a)) 

Jo 


II-  i 


2  \  BS 


«p(— Af  +  i§A/(l  —  A  —  aA))  /  1  — A 


(1  —  A)M  \l-X-a\/ 

The  value  of  A  which  gives  the  closest  bound  is  then  obtained  by  solving 


dln(CB) 


=  —  f  + 


M 

+  t - r  + 


La 


dX  '  '  (1-A-aA)2  1  —  A  (1  -  A)(l  -  A  -  aA) 

for  A.  This  result  can  be  rearranged  to  form  a  cubic  in  A 

t  A3  +  u\2  +  vX  +  w  =  0 

with  the  coefficients 

t  =  f(l  +  a)2, 

u  =  (1  +  a)  ((1  +  a)M  —  (3  +  a)r) , 
v  =  f  (3  +  2a)  —  (1  +  a)  (2  M  +  aL)  —  S, 
w  =  aL  +  M  +  S  —  f . 


=  0 


|  dp 

(C4) 

(C5) 

(C6) 


(C7) 


Equation  C6  can  then  be  solved  using  radicals.  Dividing  through  by  t  and  making  the  substitution 
A  =  rj  —  u/(3t)  gives  the  equation 


with  a  discriminant 

If  d  >  0,  there  is  one  real  root 


+  pr)  +  q  =  0  where 
o 

V  U 

p  = - „  and 

t  ,  3t2 

2 u6  uv  w 
Q  ~  27f3  “  3*2  +  J' 


d=(q/2)2  +  {p/3f. 


(C8) 


(C9) 


Ao  = 

If  d  =  0,  there  are  three  real  roots  given  by 

f  4>  +  2tt 

Aq  =  2 r  cos 


{-q  +  Vd)  +  f/i  (-9  -  Vd)  -  u/(3t).  (CIO) 


u 
3 1 


where  r  =  \J — p/3  and  </>  =  cos  1  ^3)  • 


(Cll) 
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If  w  <  0,  or  f  >  aL  +  M  +  S,  then  Ao  >  0  and  the  result  is  an  upper  bound  for  the  probability 
of  detection  near  zero.  Alternatively  if  w  >  0,  or  f  <  aL  +  M  +  S,  then  Ao  <  0  and  the  bound 
is  a  lower  bound  for  the  probability  of  detection  near  one. 


Fluctuating  target  in  correlated  clutter 


The  Chernoff  bound  for  a  fluctuating  target  is  obtained  by  integrating  the  bound  for  a  fixed  target 
in  Equation  C4  with  the  fluctuation  P(s\S,  K ), 


Cb(M,  L,  a,  S ,  K,f,  A)  =  /  CB{M,  L,a,s,f,  X)P(s\S,  K)ds 


exp  (-Af  +  ,_rlX_aX) 


0  —  \f 


(1- A)M 
1  —  A 


1  -  A 


L  zK-l 


(1- A)M  \1- X-aXj  T(K)\S 


Ii 

1-  X-aXj  T(K)\§ 
1  (IOk 


K 


-  e 


-Ks 


/ ~sd~s 


x  / 

0  — At 


sX 


s  exp  (  —Ks/ S  +  - — — — 


aA 


1  -  A 


(1  —  A)M  \1  —  A  —  aA 
.-Af  /  i  -  X 


(1- A)M  Vl-  A-aA 


K 

S 

/ 


K 


s 


ds 

X 


K 

S  1  —  A  —  aA 
A 


-K 


-K 


Kl-X-aX 


(Cl  2) 


The  value  of  A  which  gives  the  closest  bound  is  obtained  by  solving 


M 


La 


S 


dln(CB)  _  ,  _  _  _ 

dX  l-X^  {1-  X){1-  A-aA)"r(f_  A-aA)  (l  -  A  (l  +  a  +  S/K)') 

which  can  be  reamnged  to  form  a  cubic  in  A  with  the  coefficients 

t  =  f  (1  +  a)(l  +  a  +  S/K ), 

u  =  M(  1  +  a)(l  +  a  +  S/K )  —  f(  2  +  a)(l  +  a  +  S/K )  —  f(l  +  a), 
v  =  f  (3  +  2a)  -  aL(  1  +  a  +  S/K)  -  2M(1  +  a)  -  S  +  (f  -  M)S/K, 
w  =  aL  +  M  +  S  —  f. 


=  0  (C13) 


(Cl  4) 


The  solution  for  Ao  is  the  same  as  the  fixed  target  case  (Equations  CIO  and  Cl  1),  albeit  with  the 
coefficients  in  Equation  C14  replacing  those  in  Equation  Cl.  Note  that  the  coefficients  for  the 
fluctuating  target  reduce  to  those  for  the  constant  target  in  the  limit  K  — >  oo. 


Probability  of  false  alarm 

The  Chernoff  bound  for  the  probability  of  false  alarm  is 

e-Af  (  1-A  \L 

CB0(M,L,a,T,  A)  =  yi_x_aX) 


(C15) 
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which  is  obtained  by  setting  S  =  0  in  Equations  C4  or  Cl 2.  In  this  case  the  equation  for  the  value 
of  A  which  gives  the  closest  bound  is  a  quadratic  A2  +  pX  +  q  =  0  with  coefficients 


P  = 


M 

t 


a  +  2 
a  +  V 


1  —  (M  +  aL)  /f 

n  1 


(Cl  6) 


and  the  desired  solution 

Ao  =  (y~P  ~  yjp2  ~  4 qj  /2. 


(C17) 


C.2  Calculation  limits 


This  section  derives  the  lower  and  upper  limits  for  calculation  of  the  detection  probability  in 
Equation  87.  These  bounds  are  important  as  an  underflow  problem  arises  when  the  threshold,  f  is 
large  in  the  e~Tf  m/ml  terms  of  Equation  74.  As  explained  in  [Bocquet  2012],  it  is  therefore 
necessary  to  begin  the  sum  at  a  different  index  mo  >  0  so  that  the  expression  can  be  successfully 
calculated.  Shnidman  [1989]  describes  a  method  to  determine  mo  given  f  and  an  acceptable  error 
bound  e.  It  is  also  desirable  to  determine  a  priori,  a  term  mi  at  which  the  sum  will  have 
converged  to  within  e.  This  avoids  the  need  to  check  for  convergence  at  every  term  added  to  the 
sum.  These  terms  can  be  approximated  by  a  Gaussian  for  large  f ,  as  shown  in  Appendix  A  of 
[Bocquet  2012]: 


e 


m! 


exp 


(Cl  8) 


Hence  the  error  function  provides  an  approximation  to  the  sum: 


'  1  ~ra 

e-T — -  «  2  (1  +  erf(ii))  with  u  =  (mi  —  f)  /V 2mi .  (C19) 

?7l=0 


We  require 


m  i  ~m 

~i(i-erfW)<£ 

m— 0 

(C20) 

where  a  lower  bound  for  the  error  function  is  given  by  Nandagopal,  Sen  &  Rawat  [2010]: 

\/ 1  —  e-“2  <  erf(w). 

(C21) 

Thus  we  can  set 

e  =  |^l  —  \J\  —  e~u2^j  >  \  (1  —  erf(rt)) 

(C22) 

and  solve  for  u 2  to  get 

u2  =  —  In  (4£(1  —£))«—  In  (4e)  =  G. 

(C23) 

Finally,  mi  is  found  from  Equations  C19  and  C23  with  the  ceiling  function  applied  to  get  an 
integer: 


mi 


f  +  G  +  y/G{G  +  2f ) 


(C24) 
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Unfortunately,  the  approximation  in  Equation  C18  is  poor  for  m  <  f ,  particularly  if  f  is  small,  so 
the  starting  term  mo  can  be  found  using  a  different  approximation  described  by  Shnidman  [1989]: 


mo 

£■ 

m= 0 


m\ 


\  (1  +  erf (u))  with  \plu  =  \j (m0  -  |)  -  -  (m0  -  |). 


(C25) 


Here  we  require  an  upper  bound  for  the  error  function,  such  as  the  one  given  in  [Nandagopal,  Sen 
&  Rawat  2010]: 


sj i  _  e-4“2A  >  erf (n) 


so  we  can  set 


£  =  \  ^1  +  y/l  —  e~4u2rrSj  >  \  (1  +  erf(n)) . 

From  Equations  C25  and  C27, 


v2  =  =  G'  =  f  -  y2f  -  (m0  -  5)  ym,0  - 


Dropping  the  1/2  then  gives  the  solution 


m0  = 


0, 


-  y/G'(2f  -  G')  ,  f  >G' 
f  <  G'. 


(C26) 

(C27) 


(C28) 


(C29) 


The  parameter  e  is  set  so  it  is  much  smaller  than  the  precision  of  the  computer  arithmetic  (typically 
~  10-16  for  double  precision).  For  this  work,  we  use  e  =  10-30  for  which  G  ~  67.7.  Although 
we  have  from  Equation  C28,  G'  =  ttG'/4  ~  0.785G,  we  can  set  G'  =  G  without  affecting  the 
accuracy  of  the  result  provided  e  is  much  smaller  than  the  computing  precision. 
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